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Preface 


The idea for this book grew out of proposals at the APL86 con- 
ference in Manchester which led to the initiation of the I-APL 
(International APL) project, and through it to the availability of 
an interpreter which would bring the advantages of APL within 
the means of vast numbers of school children and their teachers. 


The motivation is that once schoo] teachers have glimpsed the 
possibilities, there will be a place for an “ideas” book of short 
programs which will enable useful algorithms to be brought 
rapidly into classroom use, and perhaps even to be written and 
developed in front of the class. 


A scan of the contents will show how the conciseness of APL 
makes it possible to address a huge range of topics in a small 
number of pages. There is naturally a degree of idiosyncrasy in 
the choice of topics - the selection I have made reflects algo- 
rithms which have either proved useful in real work, or which 
have caught my imagination as candidates for demonstrating the 
value of APL as a mathematical notation. Where appropriate, 
notes on the programs are intended to show the naturalness with 
which APL deals with the mathematics concerned, and to estab- 
lish that APL is not, as is often supposed, an unreadable lan- 
guage written in a bizarre character set. 
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Of course, the really important algorithms are the ones you 
program yourself, and the present collection is perhaps best con- 
sidered as a stimulus rather than an authoritative reference. 


At the time of writing, I-APL is available on the IBM PC and 
look alikes, and versions for other machines are imminent. In 
the spirit of the I-APL project, the proceeds of this book will be 
donated entirely to the furtherance of I-APL, and thereby to the 
promotion of APL as a tool for teaching. In so doing, it is good 
to recognise the tireless efforts of Anthony Camacho and Ed 
Cherlin in bringing the project to its present state of fruition, as 
well as the cooperation in J-APL publications of Linda Alvord, 
Tama Traberman and Garry Helzer, and the help and support of 
many other people throughout the world. 


Norman Thomson 

IBM (UK) Laboratories, 
Hursley Park, 
WINCHESTER. 

Hants. England. SO21 2JN 
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Chapter 1 


Introduction 


This compendium is designed to provide mathematics teachers 
with a wide range of mathematical programs of relevance to 
current syllabuses. A further aim is to provide these in small 
bulk, the combination of breadth and conciseness being achieved 
through use of APL which allows many mathematical algorithms 
to be expressed in a fraction of the number of statements that are 
required in other languages, thereby allowing the programming 
language to be truly subservient to the mathematics. 


This book is not designed to be in itself a tutorial for APL. It is 
assumed that the reader either possesses a reasonable practising 
knowledge of APL, or has studied one of the available tutorial 
introductions to the language. 


There are four essential reasons why computers should appear at 
all in the Maths classroom/laboratory, namely Computation, 
Investigation, Display, and Simulation. 


Some of the programs in this collection have the aim of 
describing an established algorithm or technique, thus fulfilling a 
computational role. Other programs are given with a view to sug- 
gesting an investigation. Here the algorithm is of no enduring 
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value of itself - its value is rather in the mathematical insight to 
be gained from writing it and observing its execution. Such pro- 
grams will be prefixed with the letter I for Investigation. Other 
programs are given which enhance mathematical insight by 
producing some numerical result in graphical form - these pro- 
grams reference functions described in Appendix 1, which are 
external to APL and depend on a supporting graphics system. 
Examples are to be found in Chapters 4 and 6. Finally, there are 
programs which carry out some specific simulation, e.g. allo- 
cating birthdays or dealing a hand at whist. The aim has been to 
range in level from school to University mathematics, and it is 
hoped that a wide range of practising mathematics teachers who 
can get their hands on a computer with APL will find something 
of relevance in these pages. 


While APL ts not unique in providing these capabilities, it has 
the overwhelming advantage of being much closer to conven- 
tional mathematical notation in its constructs than any other pro- 
gramming language available on microcomputers. This 
substantially reduces the danger of the computing tail wagging 
the mathematical dog, and of what starts out to be an investi- 
gation in mathematics becoming a wrestle with the intricacies of 
a computer system. The computer should become an unobtrusive 
presence in the classroom, available to provide immediate support 
when the progress of instruction makes it appropriate. In this 
way APL promises to be the biggest single motivator ever to 
arrive in the mathematics classroom. 


The term “program” will be used to cover several different con- 
structs. In the first place many of the programs given are accom- 
modated on a single line as “phrases.” Sometimes phrases 
embody a condition, in which case they have 3 parts, the condi- 
tion, the “action-if-true,” and the “action-if-false.” Other 
phrases are compounds of two simple phrases separated with a 
statement connector. There are thus three types of phrases, viz. 
simple phrases, conditional phrases, and compound phrases. All 
phrases are named for possible cross-reference elsewhere in the 
collection. Such names are generally function names in the APL 
sense, although when a phrase is a simple transliteration of a 
mathematical formula, as e.g. in the case of compound interest, 
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its name is taken to be the variable on the receiving side of the 
assignment statement. When it is not possible to embody an 
algorithm in the form of a phrase, it is made the subject of a 
multi-line program. 


Many of the multi-line functions follow a simple basic structure 
whose outline is 


Pla initialisation 

[2] £1:branch on condition 
[3] iterative action 

C4jJ L141 


In the interests of conciseness, the programs are devoid of the 
sort of data checks that are characteristic of commercial pro- 
grams in which such “cosmetics” often occupy a large bulk of the 
program code. The aim here is rather to use the programming 
language in ways that highlight mathematics. 


A few conventions are used to make for uniformity of style. 


(a) To help illustrate transition from mathematics to APL, 
equations and expressions in this font will are in conventional 
notation, while the heavier font is always pure APL. 


(b) In general the result of a program is denoted by Z, and the 
left and right arguments by L and R respectively. With function 
phrases however, the assignment(s) to Z are implicit, following 
the convention of direct definition form, which is described under 
(f) below. 


(c) Where it is necessary to give explicit descriptions of the argu- 
ments L and R, these usually follow the statement of the function 
or phrase. 


(d) Temporary variables are denoted by T,U,T1,T2,...; if these 
are counting variables I,J,...are used, and where context makes it 
desirable X,Y,N,H are also used as temporary variables. 
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(ec) Labels are denoted by L1,L2,... and are sequenced in order of 
appearance. 


(f) For phrases, a formulation known as “direct definition” is 
used. This is available in JI-APL where a and w are used to 
represent left and right arguments, rather than L and R. Where 
a program comprises a conditional statement together with two 
alternative statements one or other of which is executed 
according to the truth or falsity of the condition, this is presented 
in a single line in the form : 


statement conditional statement 
name: executed if : statement - executed if 
condition false condition true 


This should be read as: 


“name” 
is defined as 


unless 


in which case 


where the underlined words are the expansions of the colons. 


(g) Sometimes, in order to express a program as a phrase, or to 
group together two short similar or connected statements, the 
symbol ‘&’ is used as a statement connector. On some APL 
systems it iS necessary to write such statements on two separate 
lines, the leftmost one first. 


(h) Another form of statement compression combines a branch 
and an assignment on the same line, e.g. 


>D1,a<... 
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This should be read as “branch to LI, having set Z to ...” (See 
also Appendix 2 — Merging Branch Statements). 


(1) Index origin is | throughout. 


(j) Some algorithms are given in a recursive form, where this is 
judged to give insight into the underlying mathematics. In such 
cases the direct definition pattern under (f) above can usually be 
interpreted : 


name: recursive : stopping : stopping 
action condition action 
As above, read the colons as “is defined as,” “unless,” and “in 
which case.” 


(k) The symbol <~ is used to denote “is equivalent to.” 


(1) The APL code used conforms to the ISO APL standard with 


the addition that O¢F is used to mean “round to the nearest 
integer.” 


1.1. Graphics 


Graphics functions enable mathematical processes to be watched 
e.g. the convergence of a series or the solution of a differential 
equation. To exploit APL fully the teacher should make the 
extra effort required to realise such graphics using whatever ancil- 
lary package is available. Whereas the APL language itself is 
consistent over many interpreters, the techniques for producing 
graphics vary from one computer system to another. A few basic 
graphical facilities are all that is needed in order to generate 
useful displays. Appendix | contains a list of these described in a 
non-implementation-specific way, while Appendix 3 contains 
functions which implement them on IBM PC’s and compatible 
machines using a workspace PGRAPH supplied with I-APL. 
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1.2. Idioms 


There are some phrases which have occurred with great frequency 
in practical use of the APL language, so much so that experi- 
enced users immediately recognize them in the context of other 
programs. Some of these have little to do with mathematics, but 
are concerned with matters such as restructuring and displaying 
data — topics which cannot be altogether ignored in using APL 
for mathematics. Appendix 2 gives a list of idioms relevant 
within the scope of the present volume. 


Chapter 2 


Arithmetic and Numbers 


2.1. Basic Programs with Integers 


The index generator 1 1s a versatile tool for generating regular 
sets of integers, for example 


First R positive integers et 
First R even integers gee aot 
First R multiples of L > Dxik 


The first R odd integers are given by: 


ODDS > 14+2xiRh 
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ODDS 5 
Lo Ok 


The technique embodied in the above program can be extended 
to any problem involving series of numbers spaced at regular 
intervals, and will be considered further under the heading 
“Arithmetic Progressions” in Chapter 3. 
As a particular case consider the problem of scaling and labelling 
the axis of a graph. If the labelled points are to be from -20 to 
25 in 9 steps of 5, think of them first as 

(-20) + 0,5,10,... 45 
which is 20+5x0,19 in APL. 
This leads to the program 

AXTSLAB RL AIH C=7 RO 1) SRLS) <05aR ES] 
R : left-hand end of scale, 


right-hand end of scale, 
number of intervals 


AXTSLAB 20 25 9 
20. 15 Oe: Ob AO: 1S. DO DS 


2.2. Square and Triangular numbers 
First R square numbers : (1f)*2 or +\ODDS R 


First R triangular numbers: +\if or 2!-if 
(i.e. as binomial coefficients) 


APL can provide insight into why the first of these pairs of pro- 
grams are equivalent. The underlying process can be generalised 
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in the following way to calculate the first R Lth. powers of inte- 
gers, le. (1A) *L. 


v Z<L POWERS R;I 

as TL & Z<ii+LxR-1 

[2] 2£1:70 IF 0=I<I-1 

[3]: Hts 2e+\( (eZ) eller) s0)/Z 
V 


A trace of this function used to calculate the first 3 4th. powers 
1S: 


4 POWERS 3 
OO aa, GP M2 eo: ES 
14 15 32 65 108 
a. 167-62 


First we observe the cumulative sum of 19 omitting every 
fourth item, then the cumulative sum of this result omitting every 
third item, and finally the cumulative sum omitting every second 
item. 


I. Which numbers are both square and triangular? A suggested 
program which goes as far as the first R triangular numbers is 


SQUTRI : (TeRx*x2)/T<+\R«ik 


I. Which integers have the property that the last 3 digits of their 
squares equal the number itself ? 


V 2<SQXXX R3L 

[1] I+«RL1] & Z«10 

[2] £1:70 IF RL2)<I<I+1 

[3] 201,2<Z,(7=1000|I*2)/TZ 
v 


R : start and end points of search 
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2.3. Multiplication and other Tables 
Table generation can demonstrate the power of APL to good 
effect as in the following program for which R 1s the size of the 


required table. 


MULTAB : Ro.xRh<«iff 


Example : 
MULTAB 6 

1 2 3 & 5- 6 
2 4 6 8 10 12 
2 6 22 iS 8 
4 8 12 16 20 24 
By Or 2S: 20-225. 3.0 
6 122 28-24 30: 36 


Tables for addition, subtraction, exponentiation and other arith- 
metic functions are obtained similarly. For more on such tables 
see Chapter 5. 

The next program uses the encode function to find the row and 


column number of the first occurrence of the value R in the table 
L, read in row-major order. 


LOOKUP : 1+(pL)T 14+(,L)if 
Example : 


(MULTAB 6)LOOKUP 12 
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2.3.1. Log Tables 
A set of standard 4-figure log tables in base 10 is given by 
LOG+6 4610@(.1x9+190)°.+.01x 14110 


Anti-logarithms are : 


ANTILOG<0410*3+(.01x 1+1100)°.+ 
.001x 14+110 


2.3.2. Trigonometry Tables 
An entire set of sin/cos/tan tables can be obtained by 
1 2.30.0 


where R is a table of values. The table of values is in turn 
described by an outer product, e.g. 


(0,190)°.+.1x0,19 


to give a typical set of values from 0 to 90 degrees by 0.1. 
Allowing for conversion to radians, the trig tables are given by 


TRIGTAB<+1 2 3°.0(0#180)x(0,190)°.+.1x0,19 


2.4. lsomorphisms 


Isomorphisms between various addition and multiplication tables 
in “clock arithmetic” can readily be investigated using the outer 
product operator. For example, the addition table for integers 
modulo 4 has the same structure as the multiplication table for | 
2 4 3 modulo 4 and for 1 3 9 7 modulo 10. This equivalence is 
demonstrated by using these tables as indices to appropriate 
vectors, VIZ: 
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'ABCD'[1+4|7T°.4+7*+0 1 2 3] 
"ABDC'([5|fo.xT<1 2 4 3] 
'A B D C'(10|fo.xT<1 39 7] 


In all cases the result is 


ABCD 
BCDA 
CDAB 
DABC 


Another isomorphism concerns multiplication modulo 8 and mul- 
tiplication modulo 12: 


VA BC: De eee Ted: T8c 5-7 
'A BoC DY AO sex Teed. |G <a Ae 


In both cases the result is 


ABCD 
BADC 
CDAB 
DCBA 


2.5. Primes and Factors 


Primes up to R 
The following phrase defines the primes up to R as those of the 


first R positive integers which have exactly two factors, namely | 
and the integer itself. 


PRIMES : (2=+#0=Ro.|R)/R<«1R 
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Here is an alternative program which expresses directly the fun- 
damental property that primes are the numbers which do not 
appear anywhere in the body of the multiplication table. 

PRIMES : (~ReRo.xR)/R<141R 
Example : 

PRIMES 30 
2. 33 oP Ah AS AT eS 2G 29 


Factors of R (including R itself) 


The following phrase selects those integers which divide a given 
integer R exactly. 


FACTORS: (0=(18)|RB)/if 

Example : 
FACTORS 12 

123 4 6 12 

Factors of R (excluding R itself) 

This phrase uses the fact that no factor of R can exceed ’2R. 
FACS > (027 | A) 771.2. 5*R 

Example : 


PACS 12 
123 4 6 
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Problems involving Factors 


I. What can be said about integers with an even number of 
factors ? Here is a program which prints such integers in the 
range | to R: 


EVENFAC: (EVENFAC R-1).,Rx10=2|oFACS R: 
AR=0 = 40 


The same recursive program structure is used in the next 
example. 


I. Which are the “perfect” numbers, i.e. numbers which are 
equal to the sum of their factors ? The perfect numbers in the 
range L to R are given by: 


PERFNOS:(L PERFNOS R-1),RxiR=+/FACS R: 
Re = 1.0 


A faster way to generate perfect numbers is to observe that 
(2*R-1)x 1+2x*R is perfect provided that both R and 1+2*R 
are prime. 


I. Investigate +/+FACTORS FR for R both perfect and non- 
perfect. 


I. Find some pairs of “amicable” numbers, i.e. numbers such 
that each is the sum of factors of the other. Sometimes with 
problems like this it is sensible to leave the computer running 
indefinitely to report “interesting” cases as it discovers them, with 
the final program stop occurring on user interrupt. Here is a 
suggested program in this style. 
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V AMICNOS;I 
[1] JI<0 
[2] 21:301 IF (4+/FACS +/FACS I)#I<I+1 
[3] +01 IF T=+/FACS I 
[4] I,t+/FACS I 
[5] 51 
Vv 


I. Are there numbers which form a 3-link chain? simply add 
another +/FACS to line 2, and find out! 


2.6. HCF and LCM 


Highest Common Factor 


A simple way to calculate hcf(L,R) is to find the largest factor of 
L which is also a factor of R: 


HCF : 1t(TeFACTORS R)/T<FACTORS L 


Another algorithm for obtaining the HCF of two numbers is 
Euclid’s algorithm, which is recursive in structure: 


HUG i BUC. (OALR) |b 2 0elLon ¢ D1 


The basis of the algorithm is that if T is the smaller the two inte- 
gers L and R, then 


hef(L,R) = hef(T, max(L,R)(mod T)) . 


This process of reduction to simpler integer pairs continues until 
one of the integers is 0, in which case the hef is the other. 
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Lowest Common Multiple 


ICM : (LxR)+L HCF R 


Examples : 

12 HUC 18 
6 

12 LCM 18 
36 


2.7. Recurring decimals 


The following program obtains as a vector the cycle of digits 
which recur when the non-terminating decimal fraction L/R (L,R 
integer) 1s expanded. 


RECUR : L RECUR1 L,R 
RECUR1 : (LTtRL2)]),U RECUR1 FR 
Rl1)J=U+RC2)|7<10xL : LT+RL2] 


Examples : 


2 RECUR 7 
285714 
1 RECUR 17 
OF Be Be Bs 2 2 ER ee Gs 


2.8. Numbers in Different Bases 


The following series of problems involve conversions between 
numbers in decimal notation and their representations in other 
number bases. 


Express a positive integer R (or a vector of positive integers R) as 
a minimum length vector of digits in number base L : 
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BASE : ((L1i+ZLe1i[l/.,R)pL)tR 


Examples : 

8 BASE 21 
295 

10 BASE 23 7 
2 0 
3.7 


Express vector R as row vectors of digits of length Z[1] in 
number base L[2]: 


EXPRESS : Q(LL1]pLL2)) TR 


Evaluate a matrix R of row vectors interpreted as numbers in 
base L: 


EHVAL : L1iQR 
Examples : 


3 8 EXPRESS 21 100 


8 EVAL 3 8 EXPRESS 21 100 
21 100 


1.e. 21 1s 25, and 100 is 144,. 


Express decimal number R as a whole number part followed by 
first L digits of binary fraction : 


BINFRAC :(1+7),(2-1)BINFRAC 2x 
ITO dre 8 b= 1S 4G 


Express result of BINFRAC as a decimal number : 
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DECFRAC : R+.x2*- 1+10R 
Examples : 


5 BINFRAC 7.7 
710311 0 


DECFRAC 71011 0 
7.6875 


2.9. Roman Numerals 


Here are procedures for dealing with another number base! 


Convert Roman Numeral to Integer 


Vv Z<ARABIC R:T 
[1] T<1000 500 100 50 10 5 1['MDCLXVI'iP] 
[2] Z<—+/Tx 1x*xT>(pT)t1ivZ 

V 


R : Character string representing a Roman 
numeral 


2. Arithmetic and Numbers 19 


Convert Integer R to Roman 


V Z<+ROMAN BR31;T3U 
ae T<'MDCLXVI' & U<0 252 5 2 5TR 
[2] Z<UL1]o'M' & I<1 
[3] 01:70 IF 7<I<I+1 
C4] +02 IF UCI]=4 
[5] Z+Z,ULITI pT Lr) 
[6]  -L1 
[7] £2:+03 IF ULI-1]=0 
[8] Z<+Z,7TCI-0 1] 
[9] +L1 
(40) 232ZeC 142). 7CE-0 2] 
Pia =So4 

V 


Examples : 


ARABIC 'MCMLXXVIIT' 
1978 

ROMAN 1978 
MCMIL XXVIII 


2.10. Encoding and Decoding 


The functions T and 1 are amazingly versatile. 11fR for 
example is equivalent to +/R, and O1Ff 1s equivalent to selection 


along the last axis of R. Here are some further applications of 
these two primitive functions. 
Separate R into its integral and fractional parts : 


SEP : 0 1TA 


R : scalar, array or vector 
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phP 4.76 150 
4 150 
0.76 0 


Money Conversion 


How many 10p., 5p., and Ip. coins are required to make up a 
sum of R pence using the minimum number of coins? 


GOINS: © “O. 2 STH 
Example : 


COINS 77 


Mensuration : Conversion of Units 


Convert a vector R representing yards, feet and inches into a 
vector representing the same length in m., cm., and mm. : 


METRIC : 0 100 10T25.4x0 3 121R 
Example : 


METRIC 3 2 6 
8°30) S22 


that is, 3 yards 2 feet 6 inches is equivalent to 3m. 50 cm. 5.2mm. 
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2.11. Problems involving Base 10 Digits 


Find the sum of the digits in a base 10 integer R, (or vector of 
base 10 integers) : 


DIGSUM : +410 BASE R 
Example : 


DIGSUM 23 199 
oS 19 


I. The following program returns | if the two arguments L and R 
contain exactly the same digits (possibly in a different order). 


SAMEDIG:A/(TeU),(U+,10 BASE R)e 
T<+,10 BASE L 


Use SAMEDIG to find all the pairs of 2-digit integers such as 
(21,87) whose product 1827 consists of the same 4 digits. 


Vv RPT ;I3J 

[1] T+10 

[2] 21:70 IF 31<I<I+1 

[3] J+ |1000+L 

C4] £2:901 IF 99<J+J+1 

[5] 3+L2 IF ~(1,7)SAMEDIG IxJ 


[6] Diath 
[7] +12 
V 


I. Casting out Nines 


The following program describes a number game and uses the 
function ASK in Appendix 2. Type CAST9, and follow the 
instructions - 
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V CAST9;T 
[1] T+ASK'GIVE A NUMBER WITH ITS DIGITS 
ALL DIFFERENT: ' 
EZ] T<T,ASK'..AND ANOTHER USING THE SAME 
DIGITS =" 
[3] 'DIGIT SUM OF LARGER — SMALLER IS ', 
tDIGSUM |-/T 
Vv 


— after a while you may prefer to speed up play by using the 
simulated form CAST. 


CAST: DIGSUM |-—/10.1T7LR?R)],[1.5] 
T<- 1+R?10 


R : integer between 2 and 10 


2.12. Computer Arithmetic 


It is interesting to work out for one’s own computer the 
maximum precision with which a number is represented. One 
simple way to do this is to find the maximum value of N for 
which 


Z<2*N & Z+Z-1 


results in a value of Z which is not 1. Since Z*f 1s likely to be 
evaluated within the machine as *Rx@L , this test is valid only 
if you trust that the exponentiation and log routines are imple- 
mented using the full possible precision. If not, here is a 
program you can use which calculates the successive powers of 2 
(line 3) stopping whenever the difference between 2*7 and 
(2*7’)—1 is not computed as I. 
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VY PRECTST: 1:7 
[1] I<1 & T+2 
[2] £1:2L2 IF 1#7-T-1 
[3] TI<I+1 & T«2xT 


C4] >D1 
[5] 2£2:'MAX PRECISION IS ',(tI-1),' BITS' 
V 


Example (on one specific computer) : 


PRECTSL 
MAX PRECISION IS 56 BITS 


2.13. Counting Series Forwards and Backwards 


Another comparison which demonstrates that the arithmetical 
precision on any computer is finite involves taking a series which 
has a large number of very small terms in its tail, and observing 
the difference between counting the terms starting from the head 
and then from the tail. From a numerical point of view it is 
better to cumulate from the tail, since if the big terms have been 
added first, the addition of small terms, however many, may have 
no effect at all on the floating point sum. If the small terms are 
added first, however, their cumulative effect will make a contrib- 
ution. The following two programs do forward and backward 
summation for the sums of R terms of the series whose Ith. term 
is [*x-L. In general they will give different results, but you may 
have to look at a very low precision digit to find the difference! 


V Z2<L FWD R3L 
[1] I+Z<0 
[2] £1:70 IF R<I+I+1 
[3] 901,2<Z+1[Ix*-L 

V 
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VY Z2<«L BWD R3fI 
[1] I+R+1 & Z<0 
[2] 01:70 IF 0=I+I-1 
[3] +L1,2<Z+1*—-L 

V 


You should also observe which, if either, of these sums is the 
same as 


+/(i1R)*—-L 
Example (on one specific computer) : 


3 FWD 100 
1.20200740065967 

3 BWD 100 
1.20200740065968 

+/(1100)* 3 
1.20200740065968 


2.14. Complex Numbers 

On some APL systems complex numbers are available as simple 
scalars. Where this is not the case the complex number a+ib may 
be represented by the two-element vector (a,b). Addition and 


subtraction are given by APL vector addition and subtraction; 
multiplication and division are given by 


MUL : (—/LxR),+/LxOR 
DIV : LE&2 2oR,R MUL O 1 
Examples : 


21 MUL 3 4 
10 5 
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10 5 DIV 21 


3. C4 


For argument and amplitude, see section 6.6 - conversion to and 
from polar and cartesian co-ordinates. 


2.14.1. Complex Roots of Unity 
The complex Rth. roots of 1 are given by: 

CRU : &(2,R)p9(207) ,107T+((02)4R)x 1x1f 
Example : 


3éCRU 3 
~0.500 0.866 
~0.500 0.866 

1.000 0.000 


Chapter 3 


Algebra and Sets 


3.1. Some Basic Algebra 


This section begins with some programs designed to support 
some of the earliest algebraic ideas introduced in schools. These 


programs illustrate the use of the format function (#) to mix cal- 
culations and conversation, and also a prompting function ASK 
which is given in Appendix 2. 


Difference of two squares 


V DIFFSQ;T 
[a3 T<ASK'GIVE 2 NOS.:' 
[2] 'DIFF OF SQQ = ',(t7L1J)*2), 
bal (eT LD IKD i Ve seal Px2 
Gea '(X+Y)x(X-Y) = ',(6+/T),'x', 
(G5/( 7) = '4%/ Gs 
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Square of sum and sum of squares 


V SUMSQQ;T 
[1] T<ASK'GIVE 2 NOS.:' 
[2] UOQs. OF SUM. = (S$ t/t ys" e2 = * 
(+/T)*2 
[3] 'SUM OF SQQ = ',(eTL1i]*2),'+', 
(e7[2]*2),'=',o+/Tx2 
[4] ‘DIFFERENCE = 2x',(*T[1]),'x', 
($707) 3S" .62%/7 
Vv 


Examples : 


DIFFSQ 
GIVE 2 NOS.: 
B53 
DIFF OF SQQ = 16-9=7 
(X+Y)x(X-Y) = 7x1=7 


SUMSQQ 
GIVE 2 NOS.: 
4 3 
SQ. OF SUM = 7*2 = 49 
SUM OF SQQ = 16+9=25 
DIFFERENCE = 2x4x3=24 


Expansion of Brackets 


The next function concerns the expansion of a pair of terms each 
of the form (ax + by)(cx + dy). Ris a 2x2 matrix the first row 
of which is the value of the coefficients a,b and the second row is 
the coefficients c,d. The result 1s the coefficients of the expansion 
in the order x”, xy, y’. 


EXPAND : ((x#R),+/x#0 10R)(1 3 2] 
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Example : 


EXPAND 2 291 23 4 
3 2 8 


3.2. Roots of Quadratics 
This is written in the conditional direct definition form, with the 
3 sections in order corresponding to (1) real roots, (2) b?-4ac > 0 


test, (3) imaginary roots. 


QUAD: (-R[2]+1 1xT*.5)+2xR[1] 
O>7T+ 2 1 2+.xRxOR 
'TMAG. ROOTS:',s((-R(2]),(|/7)*.5)+2xR[1] 


R : vector of coefficients in descending power order 
The expression 2 1 2+.xRxOR is an APL version of b?-4ac. 


Where the roots are imaginary, the result is given as a 2-item 
vector a,b which is to be interpreted as (a+ib), (a-ib). 


Example : 

QUAD 12 1 
a a 

QUAD 111 


IMAG. ROOTS: 0.5 0.866 
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3.3. Matrix Operations 
Product of two (compatible) matrices L and R : Z+.xfh 
Inverse of square matrix R : ER 


Quotient of two (compatible) matrices Land R : CHR 


Solution of simultaneous linear equations : DEIR 
where the equations to be solved are Rx = L. If Lis a matrix 
with columns |,, L, .... ZEA is the set of solutions Rx = 1,, Rx = 
l,, etc. 
Example : 
R L 
sre! Comm | 4 6 
2 a a 2 9 
1 “2 4 3 9 
LER 
2 A 
“4 3 
3. 44 


i.e. Solution of x +y+z=4 isx = 2,y =-l,z = 3; 
2X - y- Z=2 
Ky eZ = 3 


| 


— change RHS to 6 -9 9, solution is x = -l, y = 3,z = 4. 
I. Confirm the following : 


(i) If A= BC in the matrix multiplication sense, then (HB)+.xA 
and A+.xIC are equal to C and B respectively. 
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(ii) A+.xHA and (f34)+.xA are both equal to UNIT itpA 
(See Appendix 2 for UNIT). 


3.3.1. Determinants 


The determinant of a square matrix R 1s given by the following 
algorithm which obtains the result as the product of the succes- 
sive topmost leading diagonal items which appear as the original 
matrix 1s reduced by progressively replacing it with the outer 
product of its leading row and column divided by the current 
topmost diagonal element. 


V Z2DET R3T 
[1] Z+RL131) 
[2] £1:2*+Zx1 1+h<«(141 14R)-((14R031)) 
o,x1+P[1;])#R01;1] 
[3] -+D1 IF 1<1+pR 
V 


Examples : 


DET 4 4pii6 


DET 2 209:-2 -7 3 


3.4. Polynomials 


In the examples in this section a polynomial is represented by a 
vector of its coefficients in descending power order. Coefficients 
of omitted powers must be represented by 0’s. 
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Evaluate polynomial R at value(s) L 

POLY : (8((pR),pL)pL)1R 
Example : 


I 2 POLY 221. 
6 13 


1.e. values of 2x? +x+3 atx = 1 and x=2 are 6 and 13. 


Addition and multiplication of two polynomials R and L 
PADD : (-—(pL)[oR)+((L4L) ,R)+(R4R),L 
PMULT : 1¥+/&(-1pL)OL°.xR,L4L 

Examples : 


2. 4.53: PADD 3.0 2° 1 
3.2 Ot 

213 PMULT 302 1 
63134 7 3 


16: (2x? Px-+ 3) + Gx? tax 1) = 3xt 2x? + 3x 4-4 and 


(2x? + x +3)(3x> + 2x + 1) = 6x°+3x*4+ 13x? + 4x?+ 7x43. 


Differentiation and integration of polynomial R 
PDIFF : 1+¥Rxd 1+19R 
PINTEG : (R+O10R),0 


Examples : 
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PDIFF 2 1 3 


PINTEG 2 1 3 
0.667 0.5 3 0 
Best Fitting Polynomial of given degree 


- see POLYFIT in Chapter 8 


3.5. Arithmetic and Geometric Progressions 
The formula for an arithmetic progression 1s: 
APS ALA Rix stay ad 
R : initial value (a), 


difference (d), 
no of terms (n) 


and its sumis +/AP . 
Example : 


AP 12 6 
i ve! 5 FP Os 

+/AP 12 6 
36 


The formula for a geometric progression has the same structure, 
and differs only in the arithmetical functions involved: 


GP : RC1)xRL21*0,1.RL3)]-1 
R : initial value (a), 


ratio (1), 
no of terms (n) 
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and its sumis +/GP . 
Example : 


GP 12 6 
14124 8 16 32 

+/GP 12 6 
63 


I. Confirm that the sums as given above conform to the 
standard mathematical formulae, namely 


S = 2a + (n-l).d (AP) 
S = a(r" - 1) (GP) 
(r - 1) 
3.6. Sets 


A set is conveniently modelled as a vector (numeric or character) 


with no repeated elements. The APL function p then gives the 
number of elements in the set, and membership of a set can be 


tested for by the function e. It is also convenient to use the func- 
tions WITHOUT and REMDUP which are defined in Appendix 
2, 


Union and intersection of two sets L and R are then given by 
UN : REMDUP L,R 


10,4 : L WITHOUT L WITHOUT R 


The best way in which to use the computer to illustrate the con- 
cepts of sets is to define a few simple vectors, e.g. 


U<+110 
A<+15 
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B<2+16 
C<«2Xi4 

and verify the various rules associated with sets, e.g. 


Distribution properties of union and intersection 


A UN B IX C <> (A UN B) IX (A UN C) 


A Ix B UN C <> (A IX B) UN (A IX C) 


De Moivre’s Laws 


U WITHOUT A IX B <> 

(U WITHOUT A) UN (U WITHOUT B) 
U WITHOUT A UN B <> 

(U WITHOUT A) IX (U WITHOUT B) 


Subsets and supersets 


If L is a subset of R, then the following subset and superset 
relations are true : 


SUBSET : A/LéeR 


SUPSET : R SUBSET L 


All subsets of a set 


ALLS > TT1x/T<Rp2 
Example : 

ALLS 3 
Or: “OQ. de dE. 
O110031%1 0 
1-20) A OTOL. 0 
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An alternative recursive version which brings the last column to 
the front is: 


ALT SCO AVE es gle 1 
T<ALLS R-1 : R=1 : 1 290 1 


All subsets of size L 
SOMESET : (L=+4T)/T+ALLS R 
Example : 


2 SOMESET 3 


ae) 
FOR 
OP PB 


This program can be used to find for a 3-item vector R 
(a) sums of all pairs of items > R+.x2 SOMESET 3 


(b) products of all pairs of items > Rx.*2 SOMESET 3 


3.7. Polynomial Coefficients from Roots 
The following function gives the coefficients of in descending 


power order of the polynomial whose roots are given by the 
items of vector R. 


PCOFS : ( 1* 14+10Z)xZ<((0,19R)°.=+4T) 
+.xRhx.*7<ALLS oR 


Examples : Coefficients of (x-1)* and (x-1)(x-2)(x-3). 


PCOFS 111 
1. “3: 4oe se 
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PCOMO ti - 2.3 


1 611 6 


This function demonstrates the use of inner and outer products 
to good advantage. 


First the columns of ALLS select all possible sets of roots. Next 
the x.* inner product converts each such set to a product of 
roots (note that V*O = 1 for all WV) 


The expression ((0,19f)°.=+#T7') gives a classification of these 
products with each row corresponding to the number of compo- 
nents (0,1,...). The +. inner product then sums the products of 
roots according to these groupings. 


Finally ( 1* 1+19Z) adjusts for the fact that the sums of the 
products of roots yield the polynomial coefficients with alter- 
nating signs. 


Chapter 4 


Series 


4.1. Recurrence Relations 


Terms of first order recurrence relation. 


Two programs are given, the first constructs the first L terms of 
the series u, = ky + k,u,., and the second prints all the terms 
using a recursive algorithm. 


V ZL SER1 R3ft 

[1] 2«RL,I<1] 

[2] 01:70 IF L<I<I+1 

[3] 72£1,2+Z,( 2+R)+.x1, 142 
Vv 


L : number of terms required 
R: Uo, ko, k 


TERM1 : L TERM1(LKR)1OL 
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L : Ko, k, 
R : Ug 


Example : 


8 SHR1 11 2 
i, ao LS: Ot, Oe: 227-255 
1 2 THRM1 1 


ou etc. 


Terms of second and higher order recurrence relations. 


The analogous functions SER2 and TERM2 deal with the second 
order recurrence relation u, = ky + ky,u,., + ku,» with initial 
values Up .U;. 


V Z2<L SER2 R3I 

[1] Z+(I+2)+Rk 

[2] £1:70 TF L<I<T+1 

[3] 3 701.2<Z,( 3+R)+.x1, 242 
V 


L : number of terms required 
R : Uo, u,, Ko, Kis k, 


TERM2 : L TERM2((kK 1+R) ,0+.x1,R 


Le ko, ky, k, 
R : Uo, U, 


Example: The first 8 terms of the Fibonacci series starting with 0 
1 are thus 
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8 SHR2 01011 
Ord 1. 2 3S 8. 1s 


and the terms of the Fibonacci series are generated by 
O 11 THRM2 0 1 


The generalisations of SER2 and TERM2 to higher order 
relations are: 


VY ZL SER R313;7T3U 
[1] Z<«(I<+L.5xpR)tR & T<—TI 
ia U<I+VR 
L3] 01:70 IF L<I<+I+1 
C4] 201,2<2,U+.x1,7+Z 
V 


: number of terms required 
ploy Wises. Kosky Koa 


ro 


THERM : L TERM((K 1+R),(1¥ 14R),.2+.x1,R 


L : k’s in ascending subscript order 
R : u’s in ascending subscript order 


I. Find the periodicity of the series u, = abs(u,.,) - u,... The 
following program produces the series for different start values. 


v Z<+L SERX R3I 

[1] Z<+R & [+2 

[2] £1:>0 IF L<I<I+1 

C3] 7201,2*«Z,(| 1+Z)-Z0 1+9Z] 
V 


L : number of terms required 
R: Up, Uy 
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4.2. Tests for Monotonicity 


The following phrases return 1 for a strictly decreasing sequence, 
0 otherwise. 


MONOD : A/ 1¥R>10R or MONOD : R<.<10R 


The following phrases return 1 for a strictly increasing sequence, 
0 otherwise. 


MONOI : A/ 1¥R<16R or MONOI : R<.2>10OR 


Note that the first versions of MONOD and MONOI could 
equivalently be expressed as inner products : 


MONOD : ( 1¥R)A.>14R 


MONOI : ( 1¥R)A.<14R 


If “strictly decreasing” and “strictly increasing” are relaxed to 
“non-increasing” and “non-decreasing” the > and < in the above 
should be changed to = and < respectively. 


4.3. Convergence 


The 1 function is valuable in helping to judge whether an infinite 
series 1S Or 1s not convergent. The general principle is to define 
T+125 say, and then do a + scan on an expression which 
defines the series. Here are a few simple examples with the con- 
ventional mathematical form for u, in brackets following : 


+\iP ( u, = I/n ) 
+\+2x*2 ( u, = I/n? ) 
+\ei7 ( u, = I/n! ) 
+\F7xT+1 ( u, = I/n(nt+]) ) 
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+\ Geta Tee i a= neem 4 
+\$1+@eT7 ( u, = I/inn ) 


Applying the DRAW function (see Appendix 1) with T as x-axis 
then gives a vivid realisation of the behaviour of the series. For 
example 


0 0 25 2 DRAW 
Pte CTED ot ST) 5 PX Pat 


gives (from top to bottom) plots of XI/n!, X1/n’?, and X1/n(n+ 1). 


4.4. Binomial Coefficients 


First L coefficients of expansion of (x+y). R may be integral 
or non-integral, positive or negative. 


BINCOEF: (0,12-1)!F 


Examples : 
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4 BINCOEF 4 
1464 
4 BINCOEF 2 
1°23 4 
4 BINCOEF 3 
1 $36 10 (triangular nos. with alternating signs) 


4 BINCOEF § .5 
th OS Oe 87S. 043125 


4.4.1. Pascal’s Triangle 


The coefficients of the terms in the expansion of (1+ x) for inte- 
gers up to R are given as columns of the matrix defined by 


PAS. + Ae. ASO. tA 
For coefficients of (1-x)® use 


PASN : ORo.!R<—iFf+1 


Examples : 

[PAS 4 [PASN 4 
i. eae GO.  @ 4 
Oa 2,3 O70: el 
0013 6 Or ©: 4-2 42 
0-0: 0" 4 4 OG oi “se. ee 
Oo 05 <0: 50: 1 BR 


(The [’s are present to force the result to integer, and thereby 
make the output display more concise.) 
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4.5. Successive differences of series 


The following function returns a vector of length (pR)—-Z 
resulting from taking differences L times. 


DIF : (14¥7)- 1+7+(L-1)DIF R : L=0: RB 


L : number of differences to be taken 
R : series in vector form 


If a table of successive differences is required, insert L}< before 
A @ Ea FOE ay sae 


Example : 

2 DIF (17)*2 
1 4 9 16 25 36 49 
o : Be ded 
22 2°22 


I. Using the variation suggested above, investigate the following 
successive differences : 


4 DIF |8 BINCOEF 4 


5 DIF |8 BINCOEF 5 


Difference tables are an excellent means of detecting isolated 
errors in series, as the following investigation demonstrates. 


I. Type the following sequence, and observe the effect of the arti- 
ficially introduced error. 


T+|20 BINCOEF 5 
0¢5 DIF T 
TL7I<TL7I+1 

O0¢5 DIF T 
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4.6. Fibonacci Numbers 


Fibonacci numbers satisfy the recurrence relationship u, = u,,, 
+ u,». Such series are generated by the following iterative 
program : 


V ZL FIB R3fI 
[1] «2 & Z<R 
[2] £1:70 IF L<I+I+1 
[C3] 701,2<Z,4+/ 242 


Vv 
L : number of terms (22) 
R : ug,U, 
Example : 
LO. TBO: a 


Oa ae 2 oO oc Be RS 2s 2S 
Add [2.1] +/ 2+Z to print the ratios of successive terms. 


8 FIB 0 1 


0 

.667 

~6 

e625 

O1142 3 5 8 13 


OOO OO F O&O 


A recursive function to output the terms of the series is 


FIB1 : R FIB1 R+KL 


where L and R are uy and u, respectively. Adding a stopping 
condition allows a clean finish, e.g. 
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FIB2 : R FIB2 R+ikL: £21000 : L 


I. (i) Define T+15 FIB O 1 and compare (2¥7)x 247 with 
dy 147*2. 


(ii) Define 7+30 FIB 1 1 and compare T[M] ACF TLN] 
and TLM HCF WN] for pairs of integers M and N between | and 
30. 
(iu) Define 

T<+60 FIB 11 
and investigate TLV] | TLNx1L60+N] for values of N in the range 
3 to 10. Describe the property of the Fibonacci series which this 


shows. 


(iv) What is the result of (NV FIB O 1)+.xPAS WN-1 for an 
appropriate small integer value of N ? 


4.7. Series relating to pi 


In each case R is the number of terms, and the function defines a 
series which converges to 7. 


Gregory's Series 7/4 = 1-37 +5'-... 
GREGORY: 4x-\+ODDS R 
Euler’s Series 7/4 =(2'+43-') - (234+3°)/3 + (2°+3->)/5 -... 
EULER : 4x—-\(4+4(43 2)°.*R)+R<ODDS RF 
Sum of reciprocals of squares of integers 


SREC2 : (6x+\#(1f)*2)*.5 
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Sum of reciprocals of fourth powers of integers 
SRECH : (90x+\#(18)*4)*.25 


I. Find the constants analagous to 6 and 90 which occur in the 
similarly defined functions SREC6, SREC8 , DAHELO 5 
SREC12. 


It is useful to show the convergence behaviour of these series by 


0 2.5 25 4 DRAW (125), (GREGORY 25), 
(FULER 25),(SREC2 25),[1.5]SRECu 25 


which graphs them in the following top to bottom order of their 
leftmost points — GREGORY, EULER, SREC4, SREC2. 


Further series which converge to multiples of x 
Product of (1 - 1/p?) for all primes p: 
PI1 : 6xx\1-+(PRIMES R)x*2 


Sums of reciprocal squares excluding multiples of 3: 
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PI2 : (6.75x+\+((043|R)/R<«1R)*2)*.5 
Sums of reciprocal squares of form 6n+ 1 or 6n-1: 
PI3 : (9x+\+(((6/R+1)€0 2)/R+1R)*2)*.5 
Sum of reciprocals of squares of integers with alternating signs: 
PT4 : (12x+\(Ro1 ~~ 1)x#(1f)*2)*.5 


Sum of quotients of factorials divided by products of odd 
numbers: 


PI5 = +\2x(!0,18)*x\ODDS R+i 
Examples : 


7 S&P L1.-30 
HeO00 “4.000 Se640- 323762 3.731. B2708 
32696. 3663S SH678 “s.074 

i Seri 2 15 
22090 2eI0s: 2976 82024 B.084 3.001 
3.072 3.082 3.088 3.094 

7 3¢PI3 30 
S000. 32052 324089 @.201 “S110. 23.215 
Selt@: iS2122 Seite Se2l26 

7 36PI4 10 
S404 2.8000. Sa2hS 34096. 34172 Gelig 
Sa Ta6” B.i260° Gulo2 82133 

7 36PI5 10 
2006) 22933. 2048: 32098 3.222. 33132 
3.137 3.139 93.141 3.141 
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4.8. Series fore 


(a) Rapidly convergent #1 : +\+!0,1f 


(b) Slowly convergent H2 : (1+4#1R)*1f 
Examples : 
EH1 6 
ly 2 ier 2a. Ziel Ob 22 IO) 2a I Led. 2e Fires 
H2 10 


22.25 2.3704 2.4414 2.4883 2.5216 2.5465 
2.9658 2.5812 2.5937 2.6042 


4.9. A Series for the square root of 2 


The following program gives the values of the fractions 1/1, 3/2, 
7/5, 17/12, .... in which each is derived from its predecessor m/n 
as (m+ 2n)/(m+n): 


VY 26Gb ROOT? Ail 
[1] T+0O & Z++/R 
[2] £1:90 IF L<I+I+1 
[3] R<+(1 2+.xRh),+/R 
C4] 2f1,2<Z,4+/R 


Vv 
L : number of terms 
R : two initial terms 
Example : 
8 ROOT2 1 1 


11.5 1.4 1.4167 1.4138 1.4143 1.4142 1.4142 
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4.10. Trig Series 


The successive series expansions of the first L terms for cosR and 
sinR are: 


COS : -\(R«eT)#!T*«2x 1410 

SIN : —\(R«T)+!T<ODDS L 
These can be compared with the values of 202 and 10L. 
Examples : 


C5: COS. 14201 
10.5 0.54167 0.54028 0.5403 0.5403 
(5 SIN 1),101 
1 0.83333 0.84167 0.84147 0.84147 0.84147 


4.11. Continued Fractions 


Evaluate the continued fraction corresponding to vector R 


CONFRAC: R[1]++R[2] : 2<oR : 
RL1)]++CONFRAC 1YVR 


The simplest continued fraction is 1+1/(1+1/(+1/(.. which 
converges to %2(1+J/5) = 1.618. Other interesting limits are 
given in the table below in which T +1 is the length of the vector 
generating the continued fraction. 


CONFRAC 1,792 ay. 
CONFRAC 1,Tpo1 2 J3 
CONFRAC 2,Tp4 J5 
CONFRAC 2,Tp2 4 J6 
CONFRAC 2,T7p1 4 J8 
CONFRAC 3,76 J10 


CONFRAC 2,793 6 iat 
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CONFRAC 2,Tp2 6 ae. 
CONFRAC 2,701 6 465 
CONFRAC 4,708 ae 
CONFRAC 4,794 8 J/18 
CONFRAC 4,7p1 8 J24 

In general 
CONFRAC N,Tp2xN J (1+? ) 


To observe the convergence, insert [+ before &L1] in the right- 
most conditional expression. 


Example : 


CONFRAC 1,602 
4 
~4167 
~4138 
~4143 
~4142 
~4142 


RF» F* RPO ND KN NO 


4.12. Interpolation 


If a function f has y, = f(x9) and y, = f(x,) the linearly interpo- 
lated estimate for f(x) is given by the formula 


Yo (X - X)) r Y; (K - Xo) 
(Xo - X;) (X; - Xo). 


A more general form of this is called Lagrange’s formula, which 
for a table R of (n+ 1) columns 


Yo Yi + Yn 
Xo XY eee Xn 
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uses the best fitting polynomial of degree n to interpolate y at a 
specified x-value L. 


Z+L LAGR R3T3L 
[i] Z<«I+0 
[2] 1:70 IF( 1+o0R)<I<I+1 
[3] T<R[2:(1 1toR)WITHOUT TI] 
[Sd ZeZeRlistixe/X/ A237) ) es 
[5] >D1 


WITHOUT : (~LeR)/L 


Example : Find the Sth. degree Lagrangian interpolation for the 
value x = 2/6 in the table T which gives values of sin x: 


T<(107),{£.5]7%«.1x16 


Ee 
0.09983 0.1987 0.2955 0.3894 0.4794 0.5646 
O.1 0.2 Oss O.4 O25 0.6 
(0+6)LAGR T 


O.5 


Linear interpolation between the two closest points gives 


(O76)LAGR O 4VT 
O.H4995 


Another technique for polynomial interpolation is to obtain 
vectors of linear interpolations on the polynomial interpolations 
of successively higher degree, each time offsetting the x-values 
(but not the y-values) by one further position. This technique is 
known as Neville’s algorithm. The best estimate at any stage is 
that for which the the tabular x-values are closest to and on 
either side of the x-value for interpolation. L is the polynomial 
degree and R is a table such as T with (0,x) appended as first 
column. 
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NEV : (-4(61 OFF(L-1)NEV R)xT-RL231]) 
+-/T+L OFF 14R(2:;] : L=0 : ivRkL1;3] 


OFF : (0,-L)+(LOR).,L.5]R 


Example : Find the 2nd degree polynomial interpolation for 
x = 7/6 1n the table T above. 


Pre e026 ar 
2 NEV 2d 
0.5049 0.5014 0.5001 0.5 


In practice, it 1s more useful to have an algorithm which carries 
on increasing the degree of the interpolation until two successive 
estimates agree within a specified tolerance. In POLI below, L is 
the required tolerance, R is defined as for NEV. The first line 
does a column reordering so that the first two x-values are closest 
to and on either side of the x-value for interpolation. If there are 
insufficient tabulated values to achieve the given tolerance, the 
result of POLI is a null vector. 


v Z<«L POLI R37313X 
[1] Rh+RL31+0,4|RL23;1J-1+RL2;) J 
[2] X<1+R[2; ] 
[3] T+0,002Z<1+R[1; ] 
[4] D1:>(L<|—-/(1tZ<«(-#(e1 OFF Z)xT-R[23;1]) 

t-f7T«(I<«I+1)OFF X),1+Z)/D1 
Ga Z+(04#9Z)+Z 
V 


Example : Interpolate table T at x = z/6 with tolerance .001. 


S001. POLL LA 


Chapter 5 


Formulae and Tables 


In using APL to do straightforward evaluations of arithmetic for- 
mulae, the input expression is usually so close to the algebraic 
statement of the formula that it scarcely merits formal description 
as a program. One of the considerable advantages of APL over 
other programming languages is that where only one algebraic 
quantity is involved in the formula this can often be generalized 
to a vector, or indeed a higher order array, to give multiple eval- 
uations with no change in APL input. Here are some examples : 


Area of a circle OR*2 
Volume of a sphere O4x (R*3) 43 
Conversion from Centigrade to Fahrenheit 32+1.8xC 


If conversion of a range of Centigrade degrees to Fahrenheit 
were required at 10 degree intervals, use 


32+1.8xC<+0,10x110 
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When more than one variable is involved in the formula, outer 
products can be used to obtain simultaneous evaluations. 


Examples : 
Volume of a cylinder : Formula OH#xRx2 
Table OHo .xRx2 
Surface area of a cone: Formula ORxL 
Table ORo . xD 
Compound Interest : Formula Ax(1+R)*W 
Table Ao,x(1+R)°o.*W 
Velocity at time T : Formula U+AxT7 


Table Uo,+Ao,xT 


5.1. Compound Interest 


In the following, A stands for Amount, R for Rate, and N for 
Number of years. The problem is to compute the value of A fol- 
lowing investment for N years at rate R. First assign values to 
A, N and R: 


A<10000 20000 30000 
ReO9 <10 +i2 
N<+5 10 15 20 
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7 O8¢CI<+Ao.x(1+h)o.*N 


15386 
16105 
17623 


30772 
32210 
35247 


46159 
48315 
92870 


23674 
25937 
31058 


47347 
opie owes, 
62117 


71021 
77812 
93175 


36425 
41772 
94736 


72850 
83545 
109471 


109274 
L2o347 
164207 


56044 
67275 
96463 


112088 
134550 
192926 


168132 
201825 
289389 
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The above formula can be generalised by making it into a func- 
tion with A,R,N as right argument, and a vector 
(oA), (oR), oN as left argument: 


Gl @UGELItR es 6 (L4bL2It bE TINA) ee (2031948 


7 03 3 4 CL 10000 20000 30000 .09 .10 .12 
SD “LOr 1S: 20 


15386 
16105 
17623 


30772 
32210 
35247 


46159 
48315 
592870 


23674 
25937 
31058 


47347 
91875 
62117 


71021 
77812 
93175 


36425 
41772 
54736 


72850 
83545 
1094714 


109274 
T2901) 
164207 


56044 
67275 
96463 


112088 
134550 
192926 


168132 
201825 
289389 
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5.1.1. Present Values 


To calculate present values it is necessary only to change “mul- 
tiply” in the above formula/function to “divide” : 


7 OtPV<+Aco.t(1+R)°.*N 


6499 4224 2745 1784 
6209 3855 2394 1486 
9674 3220 1827 1037 
12999 8448 5491 3569 
12418 7711 $4788 2973 
11349 6439 3654 2073 
TYAIS A2G72 “8235. 25253 
18628 11566 7182 4459 
17023. 9OSo 8481. S010 


5.2. Mortgage Repayments 


The problem is to compute repayments for loan A repaid over N 
years at rate R with interest compounded annually (values of 
A,R and N as above). The formula for the annual repayments in 
conventional notation 1s: 


A.r.l +r)" 
(l+r)"- 1] 


Judicious insertion of outer products again gives the table from 
the APL formula. 


Formula (#12)xAx(RxNxT) + 14+7+«(14+R) *N 
Table MORT<(412)xAo.x((Ro.xNx0)xT)+ 1+ 
T+(1+R)o.*N 
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7 2¢MORT 
Zdte2t A229 585: 203436). 91429 
219683 135%62 109.56. 97.88 
Zod TAT eho A223 5- 257 


HZ28«¢49 259.70 206.76 182.58 
439.660 27/2424 219212 195.77 
462.35 294.97 244.71 223.13 


O42573: BOF. 20: S10 et 6 2734607 
659.49 406.86 328.68 293.65 
693.52 442.46 367.06 334.70 


Thus, for example, it costs 91.29 per month to repay 10000 at 
9% over a period of 20 years. Notice how the “multiply” outer 


product balances array sizes rather than doing any direct compu- 
tation. 


5.3. Triangle Formulae 


Area of Triangle given sides 
The formula ./s(s-a)(s-b)(s-c) is given in APL as 
TARBA “2 “UX/ (y5%4/A) 80) *.5 
R : vector of lengths of 3 sides 
Example : 


TARHA 3 4 5 
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Sin Formula 


If R is a 3-1tem vector representing the angles A and B in 
degrees, and the length of side b opposite B, then the length of 
side a is given by 


SIDEA : R[3]x+/102+R 
Example : 


SIDEA 60 60 3 
3 


If R is a 3-item vector representing the sides a and b and the 
angle B opposite b in degrees, then one of the possible values for 
angle A (in degrees) is 


ANGLEA : (180401)x 10(100R[3]+180) 
x+/2+R 


Example : 


ANGLEA 3 5 90 
Clem ow | 


If R defines an impossible triangle, e.g. if a.sinB > b, then a 
DOMAIN ERROR will be reported for attempting the arcsine of 
a value greater than 1. 


Cos Formula 


If R is a 3-item vector representing the lengths of 2 sides and the 
angle in degrees between them, then the length of the third side is 
given by 


SIDE3 : ((T,200R[3]#180)+.~x 
T,x/ 2,7+2t+Rh)*.5 
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Example : 


SIDE3 3 4 90 


5.4. Longest and Shortest Journeys 


Given a distance table M for 4 towns, e.g. 


NO mM £F © 
Mb wor 
Oi 2? GO) 
OwoOFRN RS? 


the longest 2-leg journeys which can be made between pairs of 
towns is 


M[.+M 
12 9 8 15 
9 8 10 12 


A LO! do. 3g 
12° 39... “6215 


To find the shortest 2-leg journeys it is necessary to fill the 
leading diagonal of M with large numbers, e.g. 100 would do in 
this case: 


M 
100 4 6 2 
4 100 3 1 
6 3 100 g 
2 1 6 100 
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ML.+M 
4 3 7 5 
3 2 7 =«6 
710 6 4 
5 6 4 2 


For journeys of more than 2 legs the inner product is repeated 
the appropriate number of times by using “execute” to generate 


the appropriate number of [’s or L’s. 
LONGEST: #'R',(4xL—-1)o'![.+R' 
SHRTEST: #'R',(4xL-1)p'L.+R' 

L : no. of legs in journey 


R_: distance table with large value in leading 
diagonal in case of SHRTEST 


5.5. Pythagoras’s Theorem and Norms 
The formula ,/(a* + b? + c? + ... ) is given by: 


PY: (+/R*2)*.5 
R : vector a,b,c,... 


Examples: 


PY 3.4 


PY 3°45 
TeO74 


If R is generalised to be a matrix, each column of which is a 
vector of co-ordinates in n-dimensional space, an “outer product” 
of Euclidean norms is given by: 


5. Formulae and Tables 63 


V Z<NORM F;T 
[1] Z+1 107<+(QR)+.xR 
[2] Z<+((Z°.4+Z)—-2x7)*.5 


Vv 
Example : 
Ji 
O2-3-"12 
O4 5 
NORM T 
0: 5 alee: 
5 O 9.055 


1: Die0oo: 10 


Another norm is sometimes called the “taxi-cab metric” or 


L-norm. It is defined as |x, - x | + ly, - yol for points (Xo,yo). 
(x,,y,). The corresponding outer product form is: 


DNORM : ((QR)+./R)—-(QR)+. LR 
Example : 


ILNORM T 
O 7 17 
7 0 10 
17 10 0O 


5.6. Pythagorean Triples 


Pythagorean triples are sets of integers which satisfy the equation 
a* = b? + c?. The following function returns a single triple from 
scalar parameters L and R. These should be coprime to give a 
fundamental triple, i.e. one with no common factor amongst its 
terms. 


TRIPLE : TLAT+(2xLxR),|—-/T],+/T+(L.R) x2 
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The next function generates families of triples — integer R is the 
stopping criterion; the number of triples is about (but not greater 
than) ’%R(R+1). For function EUC see section 2.6. 


V Z2<PYTH F3T:3U 


[1] Tet 
[2] 01:70 IF (14+2xR)<T<+T+2 
[3] U+ 1 


[4] £2:201 IF T<U<U+2 

[5] >DL2 IF 124U FUC T-U 

[6] +02,0kU TRIPLE T-U 
V 


Example : 


Pye 4 

3 4 95 

S.-i. 

o. oZ% AS 

12 Oo od 

7 24 25 

20° 2A 29 

16 63 65 

9 40 41 

28 45 53 


Chapter 6 


Geometry and Pattern 


Most of the programs in this chapter compute arrays of numer- 
ical values representing configurations of lines and points, which 
are supplied as right arguments to the graphics programs 
described in Appendix | to produce geometrical displays on the 
screen. The functions which are used in this chapter are 
summarised here with diagrams to provide a visual reminder of 
the parameters. 


(a) DRAW - window DRAW (NxP) matrix 


N 
observns. 


and similarly for PLOT except that points are not joined. 


(b) JOIN : window JOIN (x,y) 
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This is a “draw-as-you-go” function, which is used to interleave 
calculating and plotting, e.g. when there is not enough workspace 
to store an entire array for DRAW. 


(c) LINES : window LINES (Nx4) matrix 
N Xif Yij Xo} Yao 
lines. 
(d) SKETCH : window SKETCH (Nx3) matrix 
N Oor}| x} y 
observns. 1 


Column 1 = 0 for pen-up, | for pen-down. 


(e) CIRCLES - window CIRCLES (Nx3) matrix 


N X yir 
circles 
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6.1. Parametric Plotting 


The points of a circle of radius L are defined by the single 
parameter @ where 


Lcos@ 
Lsin@ . 


X 
Ms 


I 


This leads to the following function whose result is a 2-column 
matrix whose rows are coordinates of points on the circumference 
of a circle. 


CLAC  } QL*%2 1° ,OORF180 


L : radius 
R : vector of angles in degrees 


Further examples of parametric plotting follow. 


6.1.1. Conic Sections 


Matrix of points on ellipse and hyperbola 
ELLIPSE: ((9Z)oL)xZ<1 CIRC R 


L : 2-item vector — major, minor semi-axes 
R : vector of angles in degrees 


HYPERB : REFLX((0Z)pL)*x 
Z<6 5°.00R+180 


L : 2-item vector — major, minor semi-axes 
R : vector of angles in degrees for 
one half of curve 
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REFLX : (0R),[1]Rx(pR)p1 1 
The auxiliary function REFLX returns the result of adding a 


reflection in the x-axis together with a reordering of the points 
for a continuous plot. To draw the complete hyperbola use 


“14°97 14 7 DRAW 2 1 HYPERB 5x130 


for the right hand half, then 


“14 07 14 7 DRAW 2 1 HYPERB 5x130 
for the left hand half. 
Matrix of points on parabola 

PARAB : REFLX Lx(Rx2),11.5]2xR 


L : parameter a in x=at’, y=2at 
R : vector of values of parameter t 


Example : The parabola y* = 8x is drawn by 


O 50 300 50 DRAW 2 PARAB 0,112 


6.1.2. Hypocycloids and Epicycloids 


These curves give widely varying, and in many cases aesthetically 
pleasing symmetric patterns constructed by joining successive 
points whose co-ordinates are given by the parametric equations: 
xX = cost + a.cos bt 
y = sint + a.sin bt 


V Z<«L HHPARM R 
[1] L+0(2x0,10)#L<|LxR2)-1 
[2] Z+((20L)+R[1]x20R(2]xZ),[1.5] 
(10L)+RC1)Jx10R(C2)]xZ 
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L : (number of arcs)/abs(b-1) 
R : parameters a and b 


Example : 


“4 4 4 4 DRAW 1.75 EHPARM 2 33 


The above curve connects 56 points in 7 clusters of 8. Changing 
the parameters gives considerable scope for investigation and 
experimentation. 


Hypocycloids and epicycloids are generated by a small circle 
rolling along the circumference of a larger circle — externally for 
an epicycloid, internally for a hypocycloid. The following alter- 
native function gives the parameters a more obvious physical 
meaning : 


EHCURVE: RL2] FHPARM (R(3)]+1-RL3]),1+R01) 
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R : circles per revolution, 


no. of points plotted on each circle, 
ratio of radius of revolving circle to 
that of base circle( < 1) 


RL1J>0 gives an epicycloid, R[1]<0O gives a hypocycloid. R[2] 


need not be an integer, but the product x/R[1 2] must be inte- 
gral. 


Example : 


“4 “4 4 4 DRAW FHCURVE 64 2.9375 .65 


Again, changing the parameters gives considerable variation in 
the patterns produced. 
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6.2. Envelopes 


An envelope is a curve which touches each member of a system 
of lines or curves. One of the simplest such systems is the 
“sliding ladder” curve or “astroid” : 


Vv Z+ASTROID;I 
[1] Z+O0 4pl<+.4 
[2] £41:202 IF 1<I«I+.1 
[3] Z+Z,01JIT,0 0,(1-I*2)*.5 


[4] >D1 
[5] £2:2<4 REFLECT 1 REFLECT Z 
V 


(see Appendix 1 for REFLECT) 


Lines 1-4 describe the lines in the first quadrant using Pythagoras 
theorem (line 3) applied to a ladder of unit length. Line 5 
produces the reflections in the other three quadrants. The picture 
is produced by e.g. 


“4.5 1.5 1.5 1.5 LINES ASTROID 
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6.2.1. Conic Sections 


Ellipses and hyperbolae 


Start with a base circle (say the unit circle) and a point (focus) on 
the x-axis. Join the focus to a point on the circle and draw a line 
at right angles — this line is one element of the envelope. 


v Z<«L CONIC R;T 
ea L+1 CIRC(360+L)x0,1L 
[2] Z2+(L0C3;1J]-(£20;2)]x2+L0;2])+7+R-LE31]), 
Peo)" 22 
ka 2+Z,(L0311+(L0;21x2-L0;2])#7),01.5]2 
V 


L : no. of lines in envelope 
R : x co-ordinate of focus 
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If the absolute value of the right argument R of CONIC is less 
than | the result is an ellipse, if greater than 1 it is a half of a 
hyperbola. An ellipse with foci at (.1,0) and (-.1,0) is given by 


“Ot <8: B38: DANES 35 CONLG 24 


tI KOSS 
SS 
S SA SAME ag ee 


™ %, a 
71 eee Te : 20) 
TN eet ree = 
Ter 
f, i} my 
@) , 


pee oO 4, 
eso 

>. 
SOS 


ty SSS 
SN 


The right half of a hyperbola is given by 


3) 38: 3b NES 85. CONTE 1.8 


The second half is given by changing 1.1 to -1.1 so that the com- 
plete envelope is 


3 3 3 3 LINES 
(35° CONLE 323) .Eh1135: CONTIG: 4.4. 
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R=1 gives a set of rays through the focus. 

Avoid using “regular” values for both L and R (e.g. L=36, 
R=.5) as this is likely to give a DOMAIN ERROR for division 
by 0 in line 1. 

Parabola 

A parabolic envelope is drawn by choosing a focus on the x-axis, 


and drawing rays through it to points on the y-axis. The lines at 
right angles to these rays form the required envelope. 


PARENV : REFLX 0,(RxT), ((4-RxT) x 
T<300(1L)+3xL),[1.5]4 


L : no. of lines in each half of envelope 
R : distance of focus from y-axis 


Example : 


“41 95 5 5 LINES 10 PARENV 1 
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6.2.2. Hypocycloids and Epicycloids 


The following program generates an epicycloid if R is an integer 
greater than 1, and a hypocycloid if R is an integer less than -1. 
The number of cusps is abs(R-1). 


V ZL BKPIRYP R313;73U 
[1] +1 CIRC(360+L)x0,1L 
[2] Z<0 4pI<0 
[3] 01:70 IF L<I«I+1 
C4] 2+Z.,C1)J7TLT:) .TL(L*KU=0)4+U+L|RxT; ] 
[5] >D1 
V 


These envelopes are generated by taking L equally spaced points 
on the circumference of a unit circle and joining each point I to 


the point RxI(mod L). Then use e.g. 


“2 °2 2 2 LINES 36 EPIHYP 2 


R = 2 produces the one-cusped epicycloid shown above, other- 
wise known as the cardioid. L=36 is adequate for this. R = 3 
produces the two-cusped epicycloid, sometimes known as the 
nephoid — L=72 is suggested. Use L=144 for R=4, and 
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L=216 for R=6 and above. Higher values of L will give better 
resolution but may generate WS FULL. 


To draw hypocycloids effectively, the line elements must be 
drawn outside the base circle. The function EXTEND replaces 
each chord with extensions of length L outwards from the base 
circle in each direction. 


v Z<L EXTEND R;T 
Ea PeCa/ Reo 1) side O eae 24 
[2] (ened: Oiler Lek DIA LXee 
[3] 2+2Z,01)(A#0:3 414+2x7),R0;3 4] 
V 


L=2 for EXTEND, and L=72 for EPIHYP will, together with a 
window 6 6 6 6, give good representations of a wide range 
of hypocycloids, e.g. 

6 6 6 6 LINES 2 EXTEND 72 EPIHYP 2 gives a 
deltoid: 


The cardioid can also be produced as a special case of a curve 
called a “limacon” which in turn may be drawn as an envelope of 
circles. The element of the envelope is a circle whose centre is a 
point on the base circle, and which goes through a fixed point. 
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The fixed point may be taken along the x-axis without loss of 
generality, and so the co-ordinates to be passed to CIRCLES are 


LIMACON:R,(+/((-(CpoR) 9pL,0) + 
R<«1 CIRC(360+f)x0,1f)*2)*.5 


L : x-co-ordinate of fixed point 
R : no. of points on base circle 


The special case where L is equal to the radius of the base circle 
is the cardioid. Take the window as (-(L+2), -(L+2), L+2., 
L+2). A typical picture is 

~3 3 3 3 CIRCLES 1 LIMACON 18 (for cardioid) 
Negative values of the left argument of LIMACON put the cusp 


on the other side; a value of 0 gives a figure of eight curve with a 
“knot” in the middle. 


6.3. Transformations 


Rotation about the origin in 2 dimensions 
ROTO : R+.xQ2 1°.0(0L7180)+00 0.5 
L : anti-clockwise angle in degrees 
R : 2 element vector, or 
2 column matrix, each row of which is the 
(x,y) co-ordinate pair of a point 
Rotation about a point other than the origin in 2-D 
ROT2 : T+L(3] ROTO R-T+(pR)pLL1 2] 
L : x coord of point of rotation, 


y coord of point of rotation, 
anti-clockwise angle in degrees 
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R : 2 element vector, or 
2 column matrix, each row of which 1s the 
(x,y) co-ordinate pair of a point 


Example : Sketch the triangle (1,2),(3,4),(4,3) and its image under 
rotation of 90 degrees about (1,1). 


“4 0 4 4 SKETCH 
CO; a, dt iy 4. 1. (90: ROT? FR) had 
O- dg Ae 2d: 2: 3: Ae a 


Rotation in 3 dimensions 


v Z<«L ROT3 R 
[1] Z<R 
[2] 2032 3]<£[1] ROTO 2[32 3] 
[3] ZAC¢1. 3I~CGLE2)) ROTO Zi31 3] 
C4] Zed OieG ls) FOTO Zed 2] 

V 


L : anti-clockwise angle in degrees about Oz, 
anti-clockwise angle in degrees about Oy, 
anti-clockwise angle in degrees about Ox 

R : 3 column matrix, each row of which is the 

(x,y,z) co-ordinate triple of a point in 3D 
space 


One way to draw a projection of a cube is therefore to do a 
3-dimensional rotation using ROT3 and then perform a DRAW 
on two of the resulting columns, which is equivalent to taking a 
section by one of the three co-ordinate planes. 


If we take a unit cube with one its vertices at the origin, and the 
diagonally opposite vertex as (1,1,1), an edge trace is given by the 
matrix 
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CUBE<+ 16 3p 0 


OFPrPREFPPEROOFRREOOO Oo 
PRPROoOdOrFPRREROCOOrFRRO:)O 
OOo BE Ore Be eB So Soe eS 


and a projection on the Oyz plane is given by 


= Oo OD DRAW 
(40 50 20 ROT3 CUBE)(:2 3] 
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6.4. Perspective Drawing 


PERSPEC: (RL31)Jx7),01.5]R0;2]~x 
T<DzDL—-RE;33] 


L : perspective ratio (scalar) defined by the 
mapping (x,y,z) to (xL/(L-z), yL/(L-z)) 

R : 3 column matrix, each row of which is the 
(x,y,z) co-ordinate triple of a point in 3D 
space 


Taking CUBE as defined above, a perspective drawing is 
obtained by 


4 4 dU DRAW 
2 PERSPEC ‘40°50: 20 ROT? CUBE 
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6.5. Co-ordinate Geometry in Two Dimensions 


Argument of a vector 


The following phrase returns ,/(1?7+m?) where 1 and m are the 2 
items of vector R. 


ARG : (+/Rx*2)*.5 


Signed perpendicular distance of point from line Ix +my+n 


Formula is (Ix; + my, + n)//(l? + m?). If m # 0 (ie. the line is 
not parallel to the y-axis) distance is positive if the point is above 
the line, negative if below. 


DIST : ((#,1)+.xL)+(xLl2])xARG 2+L 


L : coefficients of line equation in form 
Ix + my +n = 0 
R : 2 element vector (x,y), or 
2 column matrix, each row of which 1s the 
(x,y) co-ordinate pair of a point 


If m = 0 the function is RL1]+#/ZL3 1] with distance positive 
if the point is to the right of the line, negative if to the left. 
These two phrases could usefully be combined into a single func- 
tion. 


Image of point(s) in a line Ix +my+n = 0 (I non-zero) 


IMAGE : R-2x(|L DIST R)o.~x 
2 10 30+/b[2 1] 


L : coefficients I|,m,n of equation of line 
R : 2 element vector (x,y), or 
2 column matrix, each row of which is the 
(x,y) co-ordinate pair of a point 
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If 1 = 0 then the appropriate formula is R-—2x0,L DIST R . 


Angle between two lines 


ANG : 20(R+.xL)+x/(ARG R<2+4R), 
ARG L<24L 


returns the angle between two directed lines L and R each 
defined as a 3-item vector of coefficents of a line equation in the 
form Ix + my + n =O. Lines are determined by the ratios 1, m 
and n rather than by their absolute values, and thus there are 
infinitely many representations of L and R. If the user is indif- 
ferent to which representation is used the result of ANG will be 
one of two angles which add up to z. ANG can however be 
interpreted more sensitively by observing a convention that if 
m0 then with conventional x and y axes, the positive direction 
of a line is upward or to the right if m>0, and downward or to 
the left if m<0O. If m=0 the direction is downward if 1>0. 
ANG then returns the angle between the positive directions of 
the lines L and R. 


6.6. Polar and Cartesian Coordinates 
CTOP : Q((,x1 O/R)xARG R),[.5] 
(180+01)x 30+/6R 
R : 2 element vector (x,y), or 


2 column matrix, each row of which is 
the (x,y) co-ordinate pair of a point 


Example : 
S 

he sat 

2 aw 
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1.414 45 
“23236: "26257 


PTOC : Q((pT)p,1 0/R)~x 
T<2 1°,.00(,0 1/R)+#180 


R : 2 element vector (R,6@) or 
2 column matrix, each row of which is 
the (R,@) co-ordinate pair of a point 


For all X it is true that 
X <> PTOC CTOP X <> CTOP PTOC X 


CTOP delivers its angle in an amplitude range of -90 to 90 
degrees, and consequently, the value of the argument may be 
either positive or negative. If a positive argument with the 
amplitude ranging between 0 and 180 is required, the following 
variation should be used which is based on the formula 


k-180xxR for adding or subtracting 180 depending on whether 
R is negative or positive. 


CLOPIL. s QUARG By t.5)'CCOLS*1. O/ R} 
x180xx7)+7<+(180401)x 304+/OR 
Example : 


CLOUPL “5 
12414 45 
220 LOD «4 


It remains true that for all X 


ai eee LTO COP I. x: “or CTOPI PTOC- x 
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6.7. Patterns by Plotting Large Numbers of Points 


This 1s a computer activity which has recently become popular — 
see for example an article by Barry Martin of Aston University in 
the Mathematical Gazette, Vol. 70 No. 452 June 1986 pp.140-142 
which suggests the iterative computation of equations such as 


Xn+) — Yn 7 f(x,) 
Yn+1 — A- Xp 


with ongoing plotting of points with these co-ordinates. Here is 
a program which does this 


Vo PATIOAN Fel eTex ey 
ee X<+T+L[1] & Y<L[2] & I<0 
L2] 24:70 IF LCY4)]<I<T+1 
3 X<Y—-#F 
[4] Yabo i= 


[5] T<X 
[6] ( 4tL) PIN X,Y 
E71 >D1 

V 


L: Xo, YoA, 
number of iterations, 
elements 5 to 8 = window 
R : equation f(x) in character string form 


One suggestion is 


00 .4 250 1 11 1 PATTERN 
'(xX)x(|X)*.5! 
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i, thay 4 


Next experiment by varying L; 0 0 1.9 1000 4 44 4 
produces an interesting “chained links” effect. Then experiment 
with different functions R. Iterations of as much as 40,000 have 
been suggested to achieve some patterns. 


Chapter 7 


Calculus 


7.1. Numerical Integration 


The following functions AREAS, TRAP, and SIMPSON use the 
function AXISLAB which is given in Chapter 2, and 1s now 
repeated : 


AXISLAB: R(1]+((-/RE2 1])+R03])x0,18#L3) 
R : left-hand end of scale, 


right-hand end of scale, 
number of intervals 
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7.1.1. Upper and Lower Bounds for Integration 


AREAS produces a picture of the upper and lower bounds of 
integration given by rectangles based on the subdivision of the x 
range. (see Appendix for BOXES) 


V L AREAS BsT3U3sX 
[1] X*+AXISLAB L 
[2] CELI eC P)sb0 2161/2) BOXES 


(U,14¥7),01] (U+C 14X),0,01.5]14X), 1+VTVeR 
V 


L : end-points of x-range, 
number of points of subdivision 
R : function to be integrated as character string 


Example : O 1 20 ARFAS '(X*2)x1-X' 
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7.1.2. Trapezium Rule 


v Z<L TRAP B:X 
[1] X<AXISLAB L 
[2] Z+«((-/Z[E2 1])+Z03]) 
x+/( (101 1,(203J-1)92)x#Rh)+2 


L : range of integration, number of intervals 
R : function to be integrated as character string 


First, the values of x are calculated using AXISLAB, then the 


corresponding y-values are obtained by ¢R. The estimate of the 
integral is then the product of half the interval length and the 
weighted sum of the y-values with weights 


eee an 


— the bracketed expression immediately following the +/ is the 
APL way of describing these weights. 


7.1.3. Simpson’s Rule 


v Z<L SIMPSON B;X 
[1] X<AXISLAB L 
E23) Ret (/L (02-1) )4h031) 
x+/( (1601 1,(203]-1)p4 2)x#R)+3 


L : range of integration, number of intervals 
R : function to be integrated as character string 


This procedure is a more sophisticated version of the Trapezium 
Rule in which the weights are now 


bee 2A. Doe ee Gul 
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and the weighted sum is multiplied by one-third of the interval 


length. As before the expression in brackets following +/ is an 
APL description of the weights. 


Examples : Integrate x* from 0 to 1 using 6 Trapezium/Simpson 
intervals. 


O. 2, “ORAP "ee? 
05.3907 903 

O 1 6 SIMPSON !Xx2'! 
0. 333333 


7.1.4. Adaptive Simpson’s Rule 


For practical integration, the relevant input parameter to the 
integrating routine is not how many Trapezium or Simpson inter- 
vals should be taken, but rather the degree of precision associated 
with the answer, 1.e. repeated subdivision of intervals should take 
place until this specified degree of precision is achieved. 


Sometimes the integrand may be reasonably “level” in some parts 
of its range so that a small number of intervals will suffice, 
whereas it is relatively steep or oscillatory in others where are a 
much larger number of intervals are required. 


A technique to deal with this is the adaptive Simpson rule, where 
the range of integration is divided into 2 halves, and half of the 
overall precision is associated with each. The Simpson integrals 
for each half using 2 and 4 intervals are compared - if they are 
within the required precision, the latter is accepted as the value, if 
not, a further subdivision takes place. This is a situation to 
which recursion naturally applies. Line 2 embodies this stopping 
test; if the required precision is not achieved, lines 3 and 4 are 


executed which return the sum of 2 further calls to ADSIM. 
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v Z+L ADSIM R;T 
ey Z+(Z[1 2].4) SIMPSON R 
[2] +0 IF L[3)]>|Z-(L[1 2],2) SIMPSON R 
ES Z+((L02]-(T+0.5x-/Z[2 1]),0),0.5xZL3]) 
ADSIM R 
C4] Z<+Z+((£2[1]+0.,7),0.5xZ[3]) ADSIM R 
V 


L : lower limit of integration, 
upper limit of integration, 


absolute precision 
R : function to be integrated as character string 


Example : Integrate x> from 0 to 1 to an accuracy of 5 decimal 
places. 


O 1 .000001 ADSIM 'Xx5! 
0.16667 


7.2. Root Finding 


The problem is to find to a specified precision a root of a func- 
tion of a single variable given in the form f(x)=0. 


7.2.1. Bisection Method 


In this method the limits of search are halved at every step. It 1s 
simple, but slowly convergent. 
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Vd Bin CL Hel eT exkey 
[1] T+0 
LZ] X<Z [3 4],.5x+/Z03 4] 
[3] £1:3+53 IF L[1iJ>|—-/2tY<eR 
[4] +02 IF 0#L[2]|I<I+1 
aa ‘AT TPERATION. * (el yn" xe * CE eX 1] 
[6] D2:X<+X(T7)],.5x+/X[7T<(1+04+/xY[1 3]),3] 
wa +D1 
[38] LE3:'SOLUTION = ',((4+L)exX[1]), 

' AT ITERATION ',31 


L : precision, 
iterations between print, 
lower limit of search, 
upper limit of search, 


two left arguments for ¢ 
R : function in character string form 
for which root is to be found 


Example : Find to 4 decimal places a root of x? = 2 in the 
range (1,2). 


FPi<'(X*2)-2!' 

-00001 8 12 7 4 BISECT F1 
AT ITERATION 8 X = 1.4219 
AT ITERATION 16 X = 1.4142 
SOLUTION = 1.4142 AT ITERATION 19 


7.2.2. Iteration Method 


The principle is that the function whose root(s) are required is 
expressed in the form x=g(x), and iterations x,,, = g(x,) are 
then performed. 
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V Z2<L ITEROOT R3T3X 
L1J] x<Z[2] 
[2] >-L1 IF LL1ij<|ZL2)]-T+2kR 
[3] 7+0,2<L[2] 
L4J £1:2<(201)],07)ITEROOT R 
Vv 
L : precision, 
Start value 
R : function g(x) in character string form 
for which root is to be found 


Note that convergence occurs only when the absolute value of the 
derivative of the function gis < 1. 


Example: Find the roots of x* - 6x + 1 = 0 
2schemes (a) X,.,; = (x,? + 1)/6 


/ 


g = x/3 converges for -3<x<3 


(b) Xre 1 /(6x, 7 1) 
g’ =3/,/(6x - 1) converges for x > 5/3. 


.0001 0 ITFROOT '(1+X*2)+6' 
. 16667 
apse ase 
Feely ile sc 
pA TSO 
.005 5 ITHROOT '( 1+6xX)*.5! 
P8oO2 
eo 7956 
»I073 
: LOG 
a OL 
(OL18 
L938 
,eL 96 


OO OO CO 


On en OR OR 2 OR en O12 On| 
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7.2.3. Newton-Raphson Method 


This well-known iterative method uses the fact that the tangent 
to a reasonably smooth curve is an adequate approximation to 
the curve itself over at least a small region of x. 


v L NEWRAP B:1I3;7;71;X 
Pia T+( 1+Ri's') +R 
2 Ti<(Ri's')VR 
[3] X+Z(3] 
[4] T+0 
C5] £14:2<X-(#T7)+9T1 
[6] >L3 IF L[1i]>|Xx-Z 
EFa +02 IF 0#L(02]|I<I+1 
[8] ‘AT STHRAVION “Cet ie® XS “al Bi ek 
[9] L2:X<Z 
Edo ed 
Eta) £Os*SOLUTION = *.0 (34h) 9X )« 
C AP TPR RATION. * <5 7 


L : precision, 
iterations between print, 
initial value, 
two left arguments for ¢ 

R_ : function in character string form 
for which root is to be found, 
followed by ’;’, 
followed by derivative in 

character string form. 


Example : Find to 4 decimal places a root of x? = 2 given x=2 


as a first approximation 


FQ<' (X*2)—-232xxX' 
00001 2 2 7 4 NEWRAP F2 
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AT ITERATION 2 X = 1.5000 
SOLUTION = 1.4142 AT ITERATION 3 


The following variation estimates the derivative using the first 
element of L as stepsize and consequently does not require an 
explicit expression for the derivative: 


V L NEWRAP1 R31373X 
[1] X+Z[3] 
EZ T+<0 
[3] 21:7+o#P 
Fa X<+X+L [1] 
[5] Z<«(X-LC1])—-T+((¢#R)-T) +L [1] 
[6] 223 IF L[1]>|x-Z 
[7] +[L2 IF 0#L [2] |I<«I+1 
[8] AY TTR RATION. ©. Cer io! eS: he 

(3¥Z)sX-L [1] 
[9] D2:XxX<Z 
[10]  >L1 
Lid): £3:"SOLUTION = “ CC3VL)eXxX=L11)):, 
| AT? ITERATION. "62 


L : precision, 
iterations between print, 
initial value, 
two left arguments for 4 
R : function in character string form 
for which root is to be found 


Example : Find to 4 decimal places a root of x? = 2 given x=2 
as a first approximation 


PQ2<'(X*2)-2! 
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-00001 2 2 7 4 NEWRAP1 F2 
AT ITERATION 2 X = 1.5000 
AT ITERATION 4 X 1.4142 
AT ITERATION 6 X 1.4142 
SOLUTION = 1.4142 AT ITERATION 6 


7.3. Ordinary Differential Equations 


Three methods are given, all of which share the same program 
structure, and differ only in the core statements of the loop. Line 
10 1s optional and plots the successive points of the iterative 
sequence. 


7.3.1. Euler’s Method 


V ZL BUL R31I3N383X3Y 
Be Z<«11I+<0 
[2] X<L (1) 
[3] A#A<2(C3] 
[4] Y<L C4] 
Lod N+L(-/ZC2 1])#z 
[6] 21:70 IF N<I+I+1 


[7] Z+Z4Y 
[8] Y<Y+HxoR 
[9] X<X+H 
P10} iO: 05602] ses PLN X57 
[11] 751 
V 


L : lower limit for x, 
upper limit for x, 
step-length H, 
initial value for y 

R : differential equation in 
character string form. 
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The principle of Euler’s method is to propagate values of y using 
the iterative formula y,;,, = y; + hf’(x,y) 


The execute function is used in line 8 to evaluate the gradient f’ 
at the left-hand point (x;,y,) of the current interval. 


Example : Solve the differential equation f’ = -xy of the 
Normal curve, using x = 0, y = 2/2 as start point, and pro- 
ceeding by steps of .05 to x =3. 


A+. 05X60 
VEO Ss: 605.5 (O02). 5) 80D TaN! 


7.3.2. Mid-Point Method 
Replace lines 8-9 of EUL with 


je 
Y<Y+.5xHxeR 
X<X+.5xH 
Y<T+HxoPR 
X<+X+.5xH 


and localize T. Call this function MIDPT. 


7.3.3. Trapezium Method 


Replace lines 8-9 of EUL with 


ey 

Viey 
Y<Y+HxU<o#R 
X<«X+H 
Y<+Y1+Hx.5xU+#R 
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and localize T, U and Y1. Call this function TRAPEZ. 


The aim of MIDPT and TRAPEZ is to improve on EUL. 
MIDPT does so by estimating f’ at the mid-point of the interval 
(x,,X:41), TRAPEZ by using the average values of the gradient at 
the two ends of the interval. 


I. Investigate the differences between the three methods by 
defining 


X+.05x160 
YASCO. 2° 40545 (02) 425)MIDET t=xXxy! 
Y2?4(0. 3-40 555°(O2 )*s 5) PRAPEZ “=X xy" 


and then doing the following plots : 


O .005 3 .005 PLOT X,[1.5]Y-Y1 
O .005 3 .005 PLOT xX,[1.5]Y-Y2 
0: 20008 2 ~OO0L PLOT Xx olts5 iv 


Chapter 8 


Probability and Statistics 


8.1. Discrete Probability Distributions 


8.1.1. Binomial Distribution 
The terms in the expansion of 


(L + (1-L))® (R positive integral, 0< L< 1) 
are 


BINPROB: (L*T)x((1-L) *OT7) x (7+0,1R)!R 
Example : 


.5 BINPROB 4 
006290425. Oss7s- Oe25 Us 0625 


HIST (see Appendix 1) demonstrates the above result graphically. 
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8.1.2. Poisson Distribution 


The first L+ 1 terms for mean R are e®Rk/k! for k=0O,..L, or in 
APL 


POISSON: (*-R)x(R*«T7)+!7T<0,12 


Comparison with binomial 


The difference between a binomial series with n=10, p=.1, and a 
Poisson series with mean | is shown by 


5 POISSON 1 
0.3679. °0.9679: 0.1839--0.061371 0.01533 0.003066 


6+.1 BINPROB 10 
0.3487 0.3874 0.1937 0.0574 0.01116 0.001488 
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8.1.3. Hypergeometric Distribution 


The hypergeometric distribution describes the probabilities of a 
sample of size R containing 0,1,2,... defectives, when sampling 


takes place without replacement from a population of size Z[1] 
of which ZL2] are defective. 


HYPGEOM: ( (07) !-/L)x((7+0,18)!£02))+R!Z01] 


Comparison with binomial 


Suppose that we wish to compare the probabilities of drawing 
O,1,2,... defectives in a sample of 12 from a box of 40 items of 
which 10 are defective. The binomial and hypergeometric proba- 
bilities are given by 


5¢40 10 HYPGEOM 12 
0201048 0209776. 0.242 0.30739 0222 


Sta 25 bi NPRCD. 12 
05 03168: 0.1267 0.2323" 0.2584 0.1936 


To make the distinction more explicit, observe that the first terms 
in the above two series are 


(30/40).(29/39).(28/38) .... (19/29) 
and (30/40)'? respectively, and the third terms are 


(30C 10)-CroC a) 
4o€ 12 


and (,.C,).(0.75)*°9.(0.25)? respectively. (,C, 1s the number of 
combinations of r objects chosen from n.) 
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8.2. The Birthday Problem 


The general problem is to calculate the probabilities of finding a 
pair of people with coincident birthdays in groups of 1,2,3.,... 
people. 

We generalise this problem by defining a program which returns 
the probabilities of finding a matching pair of items in groups of 


1,2,3,... L given that each item possesses independently one of a 
set of values from | to R. 


MATCHER@ 1=<\CR=0, tii) eR 
Example : 


2423 MATCHPR 365 
0.4757 0.5073 


1.e. 23 is the size of the smallest group of people for whom there 
is a greater than 50% chance of a birthday coincidence. 


Often the object of the birthday problem is to discover how big 


the group must be before the probability of finding a matching 
pair exceeds one half. Use 


30 MATCHPR 365 


Then use the following function to simulate the number of 
matching pairs in a group of L people (R as for MATCHPR). 


MATCHES: +/~(Rifk)=19R<?L0R 


Example : 


(23 MATCHES 365),23 MATCHES 365 
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8.3. Descriptive Statistics 


The standard idiom for average of a vector as given in Appendix 
2is (+/Rf)+poR. Here is a version which also produces column 
averages if R is a matrix : 


AV : (+#4R)+1+0R 
For data in vector or matrix form, the following phrase subtracts 


the mean from all values of vector R, or column means from all 
columns when R is a matrix: 


MADJ : R-(oR)oAV R 
This is then used to calculate variance and standard deviation. 
Examples : 


QXx 
> 22. © 7 es 
2o 1 18: 27) 35 HO Be i568 


AV X 
S045 
QMADS X 
2 “3 i: at 2 2 4 


“Sa 16S 1065 SS. eS ORS. eh OF 


8.3.1. Variance 
VAR : (+4(MADJ FR)*2)+1+0R 


In the case of sampling variance, i.e. where the divisor is (n-1) 
rather than n, use 


SVAR : (+4(MADJ R)*2)+ 1+1+0R 
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8.3.2. Standard deviation 


SD : (VAR R)x.5 


Examples : 


QX 
o 2-2 Hs 6 7 se a 
25 14 18 27 35 40 32 53 


VAR X 
Bb Lobes 

SVAR X 
O05 /- LO ou 7 

SD X 
2-449 11.67 


8.3.3. Partition Values 


Median and Percentiles for vector R 


The algorithm for MEDIAN must take account of the fact that if 
the number of items is odd, the value will be that of the single 
middle item following sorting, whereas if the number of items is 
even, the median value is half the sum of the two middle items. 


MEDIAN : .5x+/RL[.5x0 1+pR<RLA4RI J 


Following the function SEP in Chapter 2, the percentile algo- 
rithm uses encode (T) in line 1 to separate a mixed number 
representing the position value of the required percentile into its 
integer and fractional parts. The fractional part gives how far 
between two items in the ordered data the required item lies, and 
the resulting calculation of the percentile value itself takes place 
in line 2. 
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VZ+L PCTILE R 
[1] Z+,0 171+.01xLx 1+ 9R+RLAR] 
[2] Z<+RCZC1JJ+Z02)x-/RLZC1J+1 0] 


Vv 
L : percentile (integer) 
R : data (vector) 
Examples : 
QX 


> 2 2 & 6. F FT & 
25 14 18 27 35 40 32 53 


MEDIAN XL31] 


25 PCTILE X([;2] 
2320225 


8.3.4. Mode and Range 


Range of a vector, or column ranges of a data matrix 


RANGE : (|4R)—-LAR 
Mode of a vector R 


The mode is a vector (possibly with just one item) of the items in 
R which have the highest frequency. An outer product is used to 
match the data items against the set of unique items. Adding 
columns counts the frequencies of the unique items, and a 
reduction following comparison with the largest such sum yields 
the (not necessarily unique) mode. REMDUP is in Appendix 2. 


MODE : (([/T)=T<+/Re.=Z)/Z+REMDUP R 
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Examples : 


QX 
Ss 2 2 & 6 7 F 8 
25 14 18 27 35 40 32 53 


RANGE X 
7 39 
MODE X{;:1] 


8.4. Random numbers from Various Distributions 
R random numbers from uniform distribution in the interval (0,1): 
RND : (?T)+T<RopL/10 


L random numbers from uniform distribution in the interval 


(ALA ds). 
RUN : R(1)+(RND L)x-/R[2 1] 


L random numbers from negative exponential distribution mean 
R: 


RNE : —Rx@RND L 


L random numbers from Normal distribution, mean R[1] and 
standard deviation RL2], using the Box-Muller formula: 


RNO : RlAJ+RL2)x((2xL RNE 1)*.5)x2002xRND L 


L random Booleans (i.e. 0 or 1) where constant probability of 
drawing a 1 1s RB: 


FBO : R>RND L 
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L random variables from the set 19f whose respective probabili- 
ties are A, and also +/R <> 1: 


RSA : +/(RND L)e°.20, 14+\R 

Frequency distribution for a sample such as the above: 
RFD : +/(ipR)°.=/[RSA R 

Examples : 


RND 5 
0.5133 0.846 0.8415 0.4154 0.4679 


4 RUN 10 20 
11.78 192-72 10.33 14.96 


4 RANE 10 
29: Teale. dee 71. 1 Oa ee8 


4 BRNO 10 3 
1HeO3 12.90 20275 8.39 


10: ABO 75 
Be a (Oe 4O? cde A a, de Ok 


LO-ROA. adh 3 Ql. 222 
de 2. A ed 2 


100 RFD .4¥ .3 .1 .2 
35 33 6 26 
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8.5. Simulations 


8.5.1. Dice/Coins,etc. 

The result of R throws of a single die is given by 
DICE : ?Rp6 

Example : 


DICE 4 
i a 


In a similar way ?Rpe2 simulates R throws of a single coin. 


Multiple dice, i.e results of L sums of throws of R dice 
MDICH : +/?(L,R)p6 
Example : 


6 MDICE 2 
OB. 0, 2 or <6 


8.5.2. Buffon’s needle 


BUFFON : +/(BND R)<20R RUN 00 .5 


This program simulates R throws of a needle of length 2 on a 
striped flag with stripes of length 2. The result is the number of 
throws in which a needle cuts a stripe. The distance x of the 
centre of the needle from the nearest stripe boundary is therefore 
uniformly distributed between 0 and 1. The angle y between the 
needle and the direction of the stripes varies uniformly between 0 
and 90 degrees. The needle cuts a stripe if x < sin y. Fora 
fixed x this happens sin-'x/z/2 of the time and hence as x ranges 
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from 0 to 1, the proportion of throws which cut a stripe is on 
average the integral of sin'x/ 2/2 from 0 to 1 which = 2/x = 
.637. 


Examples : 


(BUFFON 100),(BUFFON 100),BUFFON 100 
58 64 60 


8.6. Frequency Distributions 
FREQ : +#4(1(,R-£(1])+202])°¢.=1203] 
L : lower boundary of leftmost cell, 
cell-width, 


number of cells 
R : numeric vector 


Example : 


025 FREQ 32246779 
2.2. ob. 


Note that with the above form of FREQ, values on a cell 
boundary go into the cell on the left. If the opposite convention 


is required, replace [| with MAXI : 
MAXI : (R=[R)+/R 
Often, when data is given in the form of a frequency table, it is 


useful to transform it into a single vector with repeated values. 
Do this using replicate with the frequencies as left argument, e.g. 


S22 .-274. 2.3 
A ly Ase 2 ae 


If replicate is not implemented, use the function REPLIC in 
Appendix 2. 
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8.6.1. Stem-and-leaf Plot 


Another way of summarising data in a semi-pictorial way is the 
so-called “stem-and-leaf” plot. The following program produces a 
simple form of this plot for a vector R of positive numbers. An 
estimate V is made in line 2 of a suitable number of stems. 


V STMLEAF A3I3T3;U3V3W 
[1] T<—L10@(RANGE R)+10x10@0R 
[2] V+.1x[/R<+(7 ROUND R)x10*T 
[3] U+1+(0/-7)+/10@V 
C4] I+ 1+L/L.1xR 
[5] £4:70 IF V<I<«I+1 
[6] CCCCUS4) .O)F0TT) OLX tOxe Tt). Pls 

1 O¢WLAW<+10|(T=L.1xR)/R] 


[7] >D1 
V 
Example : 


STMLEAF 4.9 2.7 3.6 5.1 2.8 2.8 3.0 
4.6 5.2 5.8 2.2 
2.0|2788 
3.0/06 
4.0|69 
520) 128 
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8.6.2. Two-way Frequency Distribution 


When data is classified according to two variables, the function 
TWOWAY produces a 2-dimensional frequency table. 


Vv Z<+L TWOWAY R;I 
ee I+O & Z«(£[133],0)p0 
L2] £1:70 IF ££2;:3]<Z«I+1 
[3] Z+Z,L01;] FREQ(RL;2)] INRANGE 
LlL231)4+202;2)]«x 1 0+7)/RE31] 


[4] D1 
Vv 
L : matrix with two rows, each of which has 
3 items defined as per L of FREQ 
R : matrix with 2 columns - each row corresponds 
to an observation, and each column to one 
of the two classifications defined by 
the rows of L. 
Example : 
QX 


3 2 2 4 6 7 7 9 
25 14 18 27 35 40 32 53 


TWOWAY X 
0 0 


OOO: Oo: 
COON O 
ONHF OO 

Oo -O OO OS 
FOoOo;9} 

O' OO: oO 
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TWOWAY uses an auxiliary function INRANGE which for a 
numeric vector L returns a binary vector of the same length with 
a 1 if the corresponding element of L is within the range defined 
by R. 

INRANGE : (+/xLo.-R)e€0 1 


Example : 


(110)INRANGE 5 7 
Or 20-0: Or 16: A oO: 3 


If MAXI was used in FREQ change 0 1 to 0 1. A simple 
extension of this function returns those elements of L which are 
within the range R : 

WITHIN : (L INRANGE B)/L 


Example : 


(110)WITHIN 5 7 


Change the 0 1 of INRANGE to 0 for strict inclusion, or to 1 
O 1if both end-points are to be included. 


8.6.3. Scatterplots 


Scatterplots in general are given by 
window PLOT Y,[1.5]X 


where Y and X are vectors of paired values. 
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Pictorial impressions of sampling from the various distributions 
can be obtained by using PLOT or PIN. You will find yourself 
forced to use the latter if the number of points plotted 1s large 
enough to force WS FULL. 


To display a one-dimensional distribution, impose a narrow 
Normal distribution on the y-axis so that there is some sepa- 
ration of the points corresponding to the distribution of interest 
on the x axis, e.g. to display a sample of 100 random Normal 
values from N(10,3) use 


1119 1 PLOT (100 RNO 10 3),[1.5] 
100° RANO O02 


Two-dimensional distributions can be shown using exactly the 


same phrase, e.g. changing the window to 1 .06 19 .06 
gives a display of a bivariate Normal distribution. Substituting 
DRAW for PLOT gives zig-zag patterns which arise from com- 
bining different types of joint distributions (Normal vs. exponen- 
tial, Normal vs. uniform, etc.) 


8.7. Regression 


Fitting polynomials of varying degree to a set of data 
POLYFIT : LL3;2JHOZL031J]°.*0,1F 
L: Ist. column = x-values, 


2nd col = y-values 
R : degree of polynomial 


Estimates of y based on above polynomial regression 
BES? = isi) POLY L POLYFIT & 


(see Chapter 3 for POLY) 
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Residuals following above polynomial regression 


RESID : LL;2]-L EST R 


Example : 


QX 
. 2 2 & & 7 7.8 
20: BAe 27. 35 GO 32 Ss 


MK PORYE LTA. 
ete: el 2 


i.e. best fitting straight line is y = 4.452x + 7.792. 


XM ESP 4 
21.42 16.88 16.88 25.96 35.04 39.58 
39.58 48.67 
X RESID 1 
32583 (22875. 14425: 2:.040 “0.04167 0.4167 
~ 7.583 4.333 


The following function plots the least squares fit (1.e. the original 
points, the regression line and the residuals) : 


v Z<I PDRAW R;I 
[1] L DRAW (11toR),(C1.5]R EST 1 
(2) T+0 
[3] £1:70 IF (1tpR)<iI<I+1 
C4] L PIN R{I;] 
[5] EL JOIN R(I31)],(R EST 1)(CT] 
[6] >D1 

Vv 


L : window 
R : data matrix, as per L of POLYFIT 
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8.8. Correlation 


Product moment correlation coefficient 


This is sometimes known as the Pearson coefficient and is given 
algebraically by 


py (x ~ U)CY . Ly) 
J2(X- My)? JZ y-by)? 


where yu, and yw, are the means of x and y respectively. 


In APL 


CCOF : (+/(MADJ R)xMADJ L)+ 
x/(CSS R),SS L)*.5 


L : numeric vector 
R : numeric vector of equal length 


where the auxiliary function SS returns the sum of squares of a 
vector R about its mean 


SS : +/(MADJ F)x*2 


8.8.1. Covariance and Correlation matrices 
COVM : (QR)+.xMADJ R 


CORM : (COVM R)+(QT)+.xT+(1, 1tpR)o 
(+4(MADJ R)*2)*.5 


R : (NxP) data matrix with rows corresponding 
to observations and columns to variables 


116 APL Programs for the Mathematics Classroom 


The items on the leading diagonal of the covariance matrix are 
the sums of squares about the mean for each of the variables. 
The remaining items are the various cross-products of the form 


X (X - Ux)(Y - Hy) 


The correlation matrix gives the product moment correlation 
coefficients between all pairs of variables in the form of a square 
symmetric matrix. 


Example : 


QX 
8 7 6 5 4 3 2 1 
Sei2. 2 Gh. UO OE EG 
25 14 18 27 35 40 32 53 


AU S24) COFr XLs3] 
0.9531 


COVM xX 

42 42 180 
“42 48 218 
“180 218 1090 


CORM X 
1 ~0.9354 0.8413 
0.9354 1 0.9351 


“OS8413 049357 4 
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8.9. Non-parametric Tests 


8.9.1. Runs Test 


The following function counts the number of runs in a sequence 


RUNS : (RL11J=RLoR])++/R410R 


R : numeric vector 


Example : 


RUNS ‘HHHHTATTTHATTTTHHHEE ' 
7 


If there are only two types of item in the sequence, the value of 
RUNS can be used to determine whether the sequence is con- 
sistent with a null hypothesis that each item is a random drawing 


of one of the two types. 


The following is a table of critical values : 


PR = 10 20 30 40 60 80 100 120 140 160 180 200 


2 6 10 14 22 31 40 49 58 68 77 86 
3 6 11 15 24 33 42 51 60 70 79 88 
95% 8 15 20 26 37 48 59 70 81 91 102 113 
] 9 15 21 27 39 50 61 72 83 93 104 115 


The 7 runs observed in the sequence of 20 characters above are 
not inconsistent with a hypothesis of randomness. 
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8.9.2. Rank Correlation 


There are two well known coefficients which test how well two 
rankings match. (The term “ranking” is used to denote “rank 


ordering” to avoid confusion with rank in the APL sense of 
ppR.) 


To be specific, let us suppose that the marks in order of 10 stu- 
dents in French and German are given by 


FREN« 75 28 32 47 83 92 69 19 45 78 
GERM<« 74 52 41 72 80 86 73 35 67 89 


so that the rankings are: 


AAFREN 
72 3° o: G10: 6 2 4-8 
AAGERM 
732589 614 4 10 


8.9.2.1. Spearman’s coefficient 


If d, is the difference in rankings for the ith. student and N is the 
total number of students, Spearman’s coefficient is defined as 


1- 62d,7/ NCN? - 1) 
In APL, if L and R are two sets of rankings, this coefficient is 
SPEAR : 1-(6x+/(L-R)*2)+(2%*3)-T+pR 


(AAFREN) SPEAR AAGERM 
O79515 


8. Probability and Statistics 119 


The following is an abridged table of critical values for a two 
sided-test : 


N= 5 6 7 8 9 10 12 14 16 18 20 


5% 1.00 .89 .79 .74 .70 .65 .59 .54 .50 .47 .45 
10% 90 .83 .71 .64 .60 .56 .50 .46 .43 .40 .38 


Thus for N=10, the probability of the Spearman coefficient 
exceeding .65 is 2x5 = 2%%. The conclusion for the marks 
data therefore is that student performances the two subjects are 
highly correlated. 


8.9.2.2. Kendall’s Coefficient 


Continuing the example of the previous section, every one of the 
io > pairs of students can be assessed as either having consistent 
rankings between subjects (scoring +1), or as having inconsistent 
rankings (scoring -1). Here again are the rankings: 


AAFREN 
72359 10614 8 
AAGE'RM 
7325 8961 4 10 


Students | and 2 are consistent in that student 2 excels student | 
in both subjects. Students 2 and 3 however are inconsistent since 
student 2 excels student 3 in French but the reverse is true for 
German. The Kendall coefficient is 


(total plus scores) - (total minus scores) 
n&2 


In APL, if L and R are two sets of rankings, this coefficient 1s 


KENDALL: +/,((UPT pR)x(xZe.-L)xxRo.—-R) 
+2!oR 
(see Appendix 2 for UPT) 
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Example : 


(AAFREN) KENDALL AAGERM 
0.8667 


The following is a table of critical values for a two sided-test : 


N= 5 6 7 8 9 10 12 14 #16 18 20 


5% 1.00 .87 .71 .64 .56 .51 .45 .41 .38 .35 .33 
10% 600° 2/3: 402 6957: <50 47 239 230 «32 «29 <27 


Like the Spearman test, this test shows high correlation between 
the performances of the students in French and German. 


More Non-parametric Statistics 


The next three tests have the general aim of testing samples as a 
whole for consistency, i.e. the null hypothesis is that the samples 
have been drawn from the same population. For the marks 
example, they derive measures for the consistency of marking 
levels across subjects rather than correlation between subject 
marks of individual students. In each case L and R are the two 
vectors to be compared by the test, 1.e. the marks themselves, as 
opposed to their rankings. 


8.9.3. Sign Test 


STEST returns a 2-item vector, the first item being the number of 
cases in which the L value exceeds the R value, and the second 
the number of cases in which it is less than the R value. It is 
assumed therefore that the items in L and R are paired, and thus 
that L and R are of equal length. 


STEST : (4+/0/7),+/-0L7+xL—-R 
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The null hypothesis is that there is no difference in the underlying 
populations from which L and R are drawn, and hence the prob- 
abilities for each I that Z[LIJ>RLZI] and ZLIJ<R(LI] are 
both one half. 


Thus to assess the result of the test compute 
SPROB : +/(1+./L STHEST #)+.5 BINPROB ol 


which is the probability of obtaining the observed value or less at 
even luck. If this 1s sufficiently low, L and R may be judged sig- 
nificantly different. Alternatively use the following table of 
2-tailed critical values for L/Z STEST R. The following table 
gives the maximum values for which the result of STEST is sig- 
nificant at the level quoted, e.g. for N=10, the probability of the 
result being 1 or 0 is less than .025. 


N= 7 8 9 10 12 14 #16 18 20 22 24 


SVQ 0 wrt rene nn een nn nn nnn nnn enen- 
2.5% 0 0 1 1 2 2 3 4 5 5 6 
5% O dt): 2 3 4 5» 5 6 7 
Example : 


FREN STEST GERM 


FREN SPROB GERM 
0.1719 


The conclusion is thus that there is no significant difference 
between mark levels in French and German. 
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8.9.4. Wilcoxon Signed Rank Test (W-test) 


For this test the non-zero absolute differences (L-R) are ranked 
(see line 1). Each such ranking is then given a positive sign if L 
is greater than R, and a negative sign if L is less than R. The 
W-statistic is the smaller of the sums of the positive and negative 
rankings and its (1-tailed) significance can be assessed from the 
table which follows. 


V Z2<«L WIEST &;T 
Ee T<«(xT)xTRANKU |T<(L-R)WITHOUT 0 
[2] Levy (AOL) tae 

Vv 
(see Appendix 2 for WITHOUT and TRANKUV) 


eb = 7 8 9 10 12 14 16 18 20 22 24 


0.5% 0 O 1 3 7 12 19 27 37 48 61 
2.5% 2 3 5 8 13 21 29 40 52 65 81 


Differences are significant if the result of the test is below the 
values given above. With the marks data the steps in the calcu- 
lation are as follows: 


FREN 

Te 26°32 A 8S 92 69: 19° 45. 78 
GERM 

74 52 41 72 80 86 73 35 67 89 
| FREN-—GERM 

124 9 25 364 16 22 114 
AA | FREN—GERM 

195 10243 7 8 6 


The 3 cases in which a French mark exceeds the corresponding 


German mark correspond to the entries 1, 2 and 4 in the last 
vector, which total 7. Hence: 


FREN WEST GERM 
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This is significant at the 2'4% level, i.e. the probability of 
obtaining a value of 7 or lower is less than .025. The W-test is 
inherently more sensitive than the sign test since it uses values 
rather than just signs. In the present example, it leads to the con- 
clusion that there is a significant difference between the mark 
levels of French and German, as opposed to the Sign Test which 
did not. 


8.9.5. Mann-Whitney Rank Sum Test (U-test) 


This test is a non-parametric equivalent of the t-test applied to a 
pair of independent (and thus unpaired) samples. L and R are 
ranked collectively (see line 1), and the U-statistic 1s based on the 
smaller of (1) the total rankings corresponding to items in L, and 
(2) the total rankings corresponding to items in R. 


V ZL UTEST BR: 

ee T<TRANKU LR 

[2] Z<+L/((+/( pL) +7) ,+/(pL)+7)—-2!1+ (pL), ef 
Vv 


For the student marks 


AAFREN ,GERM 
gt. DBs AT DOr AO SL 2615) 8: 8. oe el. eG: be 
£2: O12 


The sum of the first 10 rankings corresponding to French 1s 95, 
and of the second 10, corresponding to German, 115. The 
U-statistic is obtained by subtracting the minimum possible rank 
sum from each. Since the numbers of French and German stu- 
dents are both 10 each minimum is | + 2+ 3+... + 10 = 
55. 
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L and R need not be of equal length for this test. When they are, 


and N = oJ, significance levels are given by the following table 
which shows for example that for N=10, the probability of U 
being 16 or less is less than .00S. 


N= 5 6 7 8 9 10 12 14 16 18 20 


S 
0.5% 0 2 4 7 11 16 27 42 60 81 105 
2.5% 2 5 8 13 17 23 37 55 75 99 127 


If L and R are of unequal length, consult more detailed statistical 
tables. Running UTEST with the marks data we have 


FREN UTEST GERM 
4.0) 


This is greater than 23 and so there is no reason to reject the null 
hypothesis that similar marking standards apply to the two sub- 
jects. Comparing this with the W-test shows that as paired data, 
the marks show a significant difference between subjects. On the 
other hand, if the marks are regarded as 2 sets of independent 
samples no such evidence 1s present. 


If all the French marks are increased by 25, and all the German 
marks lowered by 25 the U-test gives 


(FREN+25) UTEST GERM-25 
14 


The internal rankings and thus the Spearman and Kendall corre- 
lation coefficients are unchanged, but UTEST now indicates sig- 
nificant discrepancy between the mark levels in the two subjects. 
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8.9.6. Goodness of Fit 


The formula for comparing a vector L representing observations 
with an equal length vector R of expected values is 


x (O-E)?/E 
or in APL 
GOFIT : +/((L-R)*2)+R 
Significance testing is then carried out using chi-squared tables 


(see below). To apply Yates’ correction only a small adjustment 
to GOFIT its necessary. 


YATES: 2 -#7'Us.0\L7R)*2) 58 


To calculate a table of expected values for given marginal distrib- 
ution vectors L and R as in the example given below use 


EXPEC : (Le.xR)++/R 
Example : 


6 14107 43 EXPEC 61 51 38 
Boao BOe: 27 oi 
Leo: tees 2049 


In using this function it is prudent to include a test that 
(+/L)<>(+/R). 


Given a two-way frequency table R, one can go straight to the 
chi-squared value by 


CHISQ : (,R)GOFIT,(+/R)EXPEC +/L1IIR 
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Example : The following results emerged from an experiment in 
which the effectiveness of a treatment was to be judged by 
counting the number of good, bad and unaffected items in two 
sets of experimental units, one of which received the treatment 
and the other not : 


Good Unaffected Bad Totals 


Treated 46 32 29 107 
Not treated 15 19 9g 43 
61 51 38 150 


The appropriate value of chi-squared is obtained by 


CHIDO 2 3046 22 29: 25.1959 
22194 


8.10. Statistical Tables 


Normal, t, chi-squared, and F tables can be obtained using the 
equations of the appropriate probability density functions and 


integrating these using the ADSIM function given in Chapter 7. It 
is useful first to establish a routine for the Beta function with 
parameters given by the 2-1tem vector R: 


BETA : (x/!R-1)#!+/ 1,R 
Note that (4) = Vaz, so B(’2,4%) = 1. 
Also rd 2) a Y2./T, BC ¥2,'/2) = JAT, 
PQ”) = VA0AJn, = B(2'2,'2) = ¥en, etc. 
8.10.1. Normal probability density function 
Function is (1//(27))exp(- 2x7) 


NPDF : (#(02)*.5)x* .5xRx2 


R : x-value or vector of x-values 
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Example : What is the probability of a random Normal variable 
lying between O and 2 standard deviations from the mean? One 
way to find this is to use Simpson’s formula to integrate the 
Normal p.d.f.: 


O 2 .0001 ADSIM 'NPDF X' 
0.4772 


Another way to calculate Normal probabilities is to use a closely 
approximating formula for the Normal integral given in “Calcu- 
lation of the Normal distribution Function,” P.A.P. Moran, 
Biometrika, 1980, 67.3, pp.675-6: 


NINT: .5+((R#3x2*.5)4++/(((0R) 00)°.+ 
(x—-7xT29)47)x10(Ro.x(2*.5)x7*+112)43) 401 


NINT 4. 4. 2 
0.1587 0.8413 0.9772 


The inverse problem of finding the Normal Rth. percentile (or 
vector of percentiles) is solved to a high degree of approximation 
by an algorithm in “Approximations for Digital Computers,” C. 
Hastings, Princeton University Press, 1955: 


NPCT : T-((To.*0 1 2)+.xC)# 
1+ ((T+(@+(1-.01xR)*2)*.5)0.%*13)4+.xD 


where C425 515517 06802853" 10010326 
D<>1.432788 0.189269 0.001308 


Example : 


NECT “95: 27 eo ‘9925 
LeO4S. 1490-26007 
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8.10.2. Student t probability density function 
Function is 1/B(’2,’2n) (/n)(1 + (x?/n))@* D/P 


t, 1s 1/2(1 + x?) 
t, is 1/2,/2(1 + %4x?)!5 
t, is 2/n./3(1 + x?/3)?, etc. 


TPDF:+(L*.5)x(BETA T)~x 
(1+ (R*«2)4ZL)*+/T<«.5x102 


L : degrees of freedom 
R : x-value or vector of x-values 


Example : What is the probability of a random t, variable lying 
between 0 and 2 standard deviations from the mean? 


Oy 2 20001 AVSIM 2° TPF x 
0.4082 


The inverse problem of finding the percentiles of t for L degrees 
of freedom is approximated by an algorithm due to K.J. Koehler 
- see Technometrics, Vol. 25 No.1, 1983. Either L or R, but not 
both, may be vectors. 


TPCT:+ .0953+( .631414+L)+ 
(.81x (-@®Rx2-R)* .5)+.076~x 
( (R<2x1-.01xR)x(C .5+(02)*.5)xDx*.5)*F0 


Examples: 


t 2O:200 BPCT. 27.5 
Peed: LEZ: 1.979 

DEPOT. OS OFS. Oo4 5 
22 SO2 Oe god Meo oo 
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8.10.3. Gamma and Chi-squared probability density 
functions 


The Gamma p.d.f. is (x*!e>*b2)/F(a). 


GPDF : (R*xU)x(*—-TxR)+(!U«14+L—-1)+ 
(T<1+L)*1tZ 


L : shape & scale parameters (a,b) respectively 
R : x-value or vector of x-values 


X? 1s a special case of Gamma with scale parameter %. The 
parameter, degrees of freedom, is twice the corresponding 
Gamma shape parameter. 


CEDOR Ss A. Sx GEOR Af 


L : degrees of freedom 
R : x-value or vector of x-values 


The first few X? curves are: 
X2, : (1/Y(2m))x ex? 
X74: ex? 
X25: (1/J(2m))x!2e8? 
X74 ; Vyxex2 


Example : What is the probability of a random chi-squared vari- 
able lying between 0 and 11 for 5 degrees of freedom? 


O: da. 3001 ADSM "3: CPDF xX" 
0.949 


Chi-squared percentiles can be obtained by first converting the 
corresponding Normal percentile using a formula called the 
Wilson-Hilferty approximation: 


CPCT: Lx (1-(2x+9xZ)- 
((WPCT R)+#3)x(24L)*.5)*3 
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L : degrees of freedom (scalar or vector) 
R : required percentile (or vector of “% tiles) 


Examples: 


1 2 CPECT 35 
3.748 5.938 

1 CPET 30°95 
23039: 34748 


Note that accuracy falls off rapidly when R is less than about 20. 


8.10.4. F probability density function 


Function 1s gal2pb/2x(a-2)/2 
BC ‘Aa, '2b)(ax + b) x(a+ b)/2 
V Z<«L FPDF R 
ee Z<(x/L*x.5xL)xRx 14+.5x14L 


[2] Z2<Z+(BETA .5xD)x(LL2)]+2L1]xR)*.5x+/L 
Vv 


L : degrees of freedom (2-1tem vector) 
R : x-value or vector of x-values 


Some F curves are: 


Fai: 1/nQ/x)& + 1) Foy: (2x41)! 
Fog: 8(x+2)3 F 4.1): 12x(4x + 1)-25 


Example : What is the probability of a random F-variable lying 
between 0 and 5 for Fis 6)? 


O 5 .001 ADSIM '5 6 FPDF X' 
03062 
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8.11. Sample Sizes 


Consider the problem of calculating the sample size necessary to 
estimate the proportion of a population possessing a given attri- 
bute subject to a predeclared confidence level that the true but 
unknown value should be within a fraction S on either side of the 
estimate. The formula is 


P(1 - P) 
(S/Z)? 


where P is an “a priori” estimate of the population proportion 
(the worst case arises, other factors being equal, when P is 0.5 — 
the closer P moves to O or I, the smaller the required sample 
size), and Z is the standard Normal percentile corresponding to 
the required confidence interval in a double sided sense. For 
example 95% and 99% confidence levels have Z-values of 1.96 
and 2.58 respectively. 


For given ranges of P, S and Z, sample size tables are obtained 
by: 


P<«.5 .1 
S«+.01 .02 .05 
Z2*+1.96 2.58 


+SSIZE«[ (Px1-P)°o.+#(S°0.#Z)*2 
9604 16641 
2401 4161 
385 666 


3458 59914 
865 1498 
139 240 
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Planes correspond to values of P, columns to confidence levels, so 
if the anticipated P could be as much as 0.5, and the sample esti- 
mate is required to be within .02 of the population value on 
either side, then sample sizes of 2401 and 4161 are required for 
95% and 99% confidence levels respectively. 


If the population size N is not indefinitely large, the sample size 
reduces according to the formula 


P(1 - P) 
(S/Z)2 + P(1-P)/N 


The convenient way to handle this is to reciprocate the previous 
table, adjust, and reciprocate: 


[+(+SSIZE)°.++N<5000 10000 


3289 
3845 


1623 
2272 


358 


988 


2045 


2726 


738 
1293 


136 
230 


49QQ 
6247 


19337 
2909 


3714 


625 


20/0 


3747 


47 
1303 


138 
235 
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An extra dimension corresponding to N has been introduced. It 
features as columns in the above display so that for example 2401 
should be adjusted to 1623 and 1937 for populations of 5000 and 
10000 respectively, and 4161 should be adjusted to 2272 and 
2939. In interpreting a table of this sort the user must assess 
carefully the validity of the assumptions which underlie the basic 
formula, particularly those concerning random sampling. 


Chapter 9 


Combinatorics 


Combinatorial algorithms deal with the arrangement of objects in 
cells. The number of such algorithms is legion due to such con- 
siderations as : are the objects ordered within cells; are the cells 
themselves ordered; are empty cells allowed; how many kinds of 
distinguishable objects are present. The whole range of such 
algorithms would fill a book in its own right — a few of the 
most basic ones are given here covering permutations (with 
derangements as a subset), combinations, selections with repe- 
tition, compositions and partitions. 


9.1. Permutations in Lexical Order 


The following algorithm gives a list of all permutations of first R 
integers in lexical order. 


The first step is to define a function which transforms an integer 
R into its representation in a number base whose elements are 


‘Oil 


FACTINT : (01+1L)TR 
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Example : 


3 FACTINT 21 
oo Aa 


1.e. 21 is represented in 3-factorial notation by the vector 3 1 1 
since 21 = (3x3!) + (1x2!) + (1xl!). Provided that the succes- 


sive elements of such a vector do not exceed 1M, then each of 
the integers 1+1!M+1 is uniquely represented. 


Suppose that we want to obtain the Ith. permutation of iW 
when the permutations are arranged in lexical ordering, i.e. the 
ordering which for N=4 begins 


1234 1243 1324 


A possible procedure is to start with iN alongside the 
(N-1)-factorial representation of Z, the successive digits of which 
are then used as indexes in origin 0 to pick out the elements of 
the permutation, always removing the chosen element of iM 
before making the next selection. 


For example with W=4, J=10 the 3-factorial representation of I 
is | 2 0, and so we start with 


1234 alongwith 120. 


Now use | from 1 2 0 to select 2 as the first element of the per- 
mutation leaving 


134 alongwith 20. 


Now use 2 from 2 0 to select 4 as the second element of the per- 
mutation leaving 


13 alongwith 0. 
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Use 0 to select 1 as the third element leaving 3 as the final 
element to give the required permutation as 2 4 1 3. 


If R is a factorial number and L is a vector to be permuted 
which may be any vector (character or numeric) of length at least 


one greater than pf, this algorithm may be expressed as 


LEXPERM: LUT], ((2LT*«1+1thIJ4Z)/LZ) 
LEXPERM 14R : O=0R : L 


Examples : 


1 2 3 4 LEXPHERM 1 2 0 
24 1 3 

'ABCD' LEXPERM 1 2 0 
BDAC 


This may now be combined with FACTINT to give a program 
which returns the Rth permutation of the vector iL. 


PERM : (1L)LEXPERM (L-1)FACTINT R-1 
so that we have for example 


4 PHARM 21 
ee 


A complete list of the first R permutations of order L is given by 


PERMS : (L PERMS R-1),l1]L PERM R 
R=0 : (0,L)90 
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3 PERMS 6 


WO WN N FP F&F 
NR RFP @®D FP W ND 
FY NM RFP WD DN W 


A permutation R can be converted into a binary form by 
BINPERM : (1ipRf)°.=R 

Examples: 
3 PERM 4 


BINPERM 3 PERM 4 


1.e. the 1’s in the columns are located in rows 2,3 and 1 respec- 
tively. 


9.2. Derangements 


A derangement is a permutation in which none of the items are 
in their natural place. Derangements of order L can therefore be 
selected from permutations by modifying PERMS to: 


DERANGE : (L DERANGE R-1),11] 
((~-v/T=1L) ,L)eT<L PERM R 
R=0 : (0,L)p0 
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The complete list of derangements of order 3, for example, is 


3 DERANGE 6 
23° a 
o a 2 


The number of derangements of a given order is 
NUMDER : (!Rf)x-/+!0,1f 
I. R people each pick up a hat from a cloakroom in which the 


lights have gone out. How many people have a wrong hat? Sim- 
ulate this situation with 


HATS : +/(.f)#R?R 


Repeat this experiment L times and find the total number of 
instances of a wrong hat : 


RPTHATS: (HATS F)+(£-1)RPTHATS R: L=0: 0 


You should observe that the average number of people with a 
wrong hat tends to R-1, e.g. 


(100--RPTAATS 6) 2100 
Bid 


Now investigate the proportion of experiments in which there is 
the hats are totally mismatched by changing (HATS R) to 
(R=HATS FR) in RPTHATS: 


(100 RPTHATS 6)+#100 
0.34 


This result should tend to l/e = 0.3679 as L and R increase. 
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9.3. Combinations 


To obtain all combinations of L integers out of R, we first revisit 
the function SOMESET in Chapter 3 which returns a matrix 
whose columns represent combinations in the form of binary 
integers. 


2 SOMESET 4 
O00 2b ake St 
Ors ae Oe Gy cat 
A 20) A. Os. Be 
a. A Or 2, 0 - 


The first column represents the combination 3 4 since the 1’s are 
in rows 3 and 4, and similarly for the other columns. 


The function BITONUM carries out this transformation 
recursively column by column: 


BTONUM: (RCE:11]/11+pR),BTONUM 0 1VR 
= 1toR :(,R)4(pR)p10,R 


Example : 


BTONUM 2 SOMESET 4 
322111 
4uU3S 4 3 2 
A function to enumerate as rows all the combinations of L 
objects out of R is thus 


COMBS : QBTONUM L SOMESET R 
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Example : 
2 COMBS 4 


RRP RP BM DY WH 
NM WO FW fF 


9.4. Selections 


An algorithm for all selections of L integers out of R, with repe- 
titions demonstrates a further application of the encode function. 
Note that the number of such selections is R'. The selections 
themselves are the columns of the result. 


SELECT : 1+(LoR)T 1+ 1R*L 


A list of all 3-letter words which can be made from the alphabet 
‘AB’ is 


'"AB'L3 SELECT 2] 
AAAABBBB 
AABBAABB 
ABABABAB 


If the selections are regarded as internally unordered use 


QREMDUPM RSORTU &2 SELECT 3 
dA de Zt GZ: 8 


gears se 
(see Appendix 2 for REMDUPM and RSORTU) 
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The number of such selections is 2!Z+R-1 . The fact that this 
is a number of combinations means that the above list of 


numbers could also have been derived from 2 SOMESET 4 (see 
previous page). Within each column interpret the 0O’s as 
“dividers” betweens 1’s, 2’s, 3’s,. Thus column | with 2 dividers 
at the start represents no 1’s, no 2’s, and two 3’s. Column 2 1s 
no l’s, one 2, one 3, and so on. The function BTONUM can be 
amended to implement this: 


BTONUM1: (1+R0311/+\~RE31]), 
BTONUM1 0 1VR : 1= 1t+oR 
1+(,R)4+\~R 


giving 


BTONUM1 2 SOMESET 4 
oa 2 de ae 
3 Oo -2 


9.5. Compositions and Partitions 


A composition of an integer L is an allocation of L indistinguish- 
able items to an unspecified number of non-empty ordered cells. 
The compositions of 4 are thus 


—- 


Bee KE HNL WwW 
poet eet ND et AQ ee 
m= KQ 


A single composition of L is given by 


COMPOSN: 1 DIF 0,((((£-1)p2) TR) ,1)/i2 
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where R is an integer in the range | to 2'', and DIF 1s as given 
in Chapter 4. A matrix of the first R compositions of L is given 
by 


COMPS : (L COMPS R-1),L1]LtL COMPOSN R 
R20. = (0-5 ):00 


To display the compositions one by one as shown above, insert 
[k before L COMPOSN R. 


If compositions such as 2 1 1 and | 2 1 which differ only in 
internal ordering are judged to be equivalent, we talk of the par- 
titions of R. Each partition corresponds to a way in which R can 
be expressed as the sum of integers less than or equal to R. Par- 
titions can be obtained from compositions by 


REMDUPM RSORTD COMPS 4 


S200 
DZ 30 
2110 
de os a oe 
400 0 


9.6. Latin Squares 


Latin squares are combinatorial patterns of letters or numbers 
which are much used in the Design of Experiments. In the 
square of order R each row and column contains exactly one 
occurrence of the integers 0,1,....R-1). The number of distinct 
Latin squares of order R grows rapidly as R increases — the 
program below gives the simplest pattern which can be obtained 
for all integers greater than 1. 


LATIN : R| (oT)e°.-T<ih 
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LATIN 4 
2 2 1.20 
yams Oa 6 
A O43; 2 
Os 22. 


9.7. Magic Squares of Odd Order 


Magic squares are patterns of numbers which belong mainly to 
the world of recreational mathematics. There are many vari- 
ations on the basic theme which is that all rows, columns and 
two major diagonals add up to the same total. As with Latin 
Squares, we give a simple algorithm which produces a magic 
square of any odd order 2R +1. 


MAGIC : Te(-oT+1R)O(L.5xR)O(20R) pi k*2 


MAGIC 5 
17 24 1 8 15 
23 5 7 14 16 
4 6 13 20 22 
POA AO 2. 8 
La Or Zo. 2 


Chapter 10 


Games and Miscellaneous 


10.1. Deal a Hand at Whist 


The following program to deal whist/bridge hands exhibits yet 
another use of encode, this time to assign a random drawing 
from the 52 cards of a pack to a suit/value pair. The assignment 
is made numerically in line 1, and card values are assigned by 
indexing in line 6 which also sorts each hand in suit order. 


V: AAND eT eet 
[1] T<14+0 1374 139 1452752 


[2] T<0 
[3] 01:70 IF 4<JI<+I+1 
C4] J<0 & '! 


[5] £2:901 IF 4<J<+J+1 
[6] 'SCHD'(CJ],'-', '234567897IQKA ' 
LULVUSCdi=P Lies) 77232 3) 
[7] +12 
V 
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Example : 


HAND 


S-QT3 
Cars 
H-AJ872 
D—K644 


p= 7 
C-A9842 
BES 
D-Q872 


S-A9862 
C=Kb 
H-K8 64 
D-AQ 


S-K54 
C-QI75 
H=O3 
D=-JT53 


10.2. Chessboard 


This program returns a chessboard pattern made up of the char- 
acters ’*’ and ‘blank’. 


CHBOARD: 'x '[1+(CHVECT B[1])eo.= 
CHVECT FLO) .) 


R : 2-item vector giving size in characters 
of one of the 64 squares of the board 
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The auxiliary function CHVECT returns alternate sequences of 
I’s and 0’s, each of length R. 


CHVECT : 2|L( 1+18xR)+R 


10.3. Mastermind 


The six colours of the standard Mastermind game are represented 
as the first six positive integers. The computer sets a pattern (line 
1), and alternately prompts and scores the user until a correct 
pattern match is achieved. A “black” is a perfect match, colour 
and position, a “white” is a match for colour only. 


V Z<MM;T3U 
[4] Z<4?6 
[2] £1:7*+ASK'YOUR GUESS?! 
[34 Si? LE AJ TSS 
C4] 'WHITES = ',(0[((+/TeZ)L+/ZeT)-U), 
' BLACKS = 1. 0US+/7=Z 


rs) «S64 

[6] 02:'WELIL DONE — CORRECT ANS. IS :-' 
Vv 

Example : 
MM 

YOUR GUESS? 

2 oO Hh 5 


WHITH =. 25. BLACKS = 1 etc. 
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10.4. Life 


The game of Life is played on a rectangular board whose initial 
dimensions are provided by the player in line 2. A set of initial 
occupied cells is determined by the player giving a series of row 
and column co-ordinates terminating with a null line (lines 3 and 
4). Then the number of generations required must be given (line 
6). 


The rules of Life are embodied in line 10, namely that at each 
new generation, each cell, whether occupied or not makes a 
count of its neighbours. A “neighbour” is an adjoining cell — 
horizontally, vertically or diagonally — so that each cell on the 
board has 8 neighbours — this is reflected in the 8 subexpressions 
in lines 8 and 9. Any cell with 3 neighbours becomes or remains 
occupied, while an occupied cell with 2 neighbours also remains 
occupied. (Occupied cells with more than 3 neighbours are 
assumed to die from overcrowding, those with less than 2 to die 
from boredom!). 


a a ee Laas ee 
[1] T<BORDER(T+ASK'SIZE OF BOARD?')p0 
[2] 'STARTING FIGURE '! 
[3] £1:9L2 IF 0=pU+ASK'NEXT CELL? ' 
C4 | See | SUE ed 
CS) 22 «eh eee 
[6] 370 IF 0=pU<ASK'NO OF GENERATIONS? ' 
[7] £3:7*L2 IF 0>U<U-1 
[8] T1<(10e7)+(10 107)+(167)+1>10e7 
[9] 71<714+('160 107)+(° 167)+( 16107)+ 107 
[10] %<(71=3)V(71=2)aAT=1 
[at)}. 363 
V 
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For BORDER see Appendix 2. The current state of the board is 
displayed as a character matrix in line 5. This game became the 
subject of an extensive literature a few years ago, see e.g. Martin 
Gardner’s articles in the Scientific American Vol 223 no. 10 (Oct. 
1970) and Vol 224 No. 2 (Feb 1971). 


10.5. Recursive Algorithms 


10.5.1. Tower of Hanoi 


At the start L rings, all of different sizes, are piled on the left- 
most of 3 pegs (peg 1) with the rings in increasing order of size 
going down the peg. The object is to repeat this configuration 
on the rightmost peg, moving only one ring at a time, and 
without ever allowing a ring to rest on a ring smaller than itself. 
The following program achieves this with the minimum number 
of moves. 


v L HANOI R 
[1] +0 IF L=0 
EQ (L-1)HANOI RL1 3 2] 
[34 Cob i OR Le "er G2 
[4] (L-1)HANOI RL3 2 1] 
V 


L : number of pegs to be moved 

R : vector of numbers representing source, 
target and spare pegs respectively 
(say 1 2 3) 


The principle is that once the first N-1 rings have been success- 
fully transferred to the spare peg (line 2), the Nth. ring is then 
transferred to the target peg (line 3), and the first N-1 rings are 
then moved on top of it (line 4). 
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10.5.2. Ackerman’s Function 


If you have not met it before, this innocent-looking function will 
amaze you by its capacity to gobble up computer time in evalu- 
ating the function for some small values of its (positive integer or 
zero) arguments. The function is defined first in algebra, and 
then in APL : 


f(0,R) = R+1 
f(L,0) = f((L-1),1) 
f(L,R) = f((L-1),f(L,R-1)) 


V Z<L ACKER R 

ea PIGD- Le OM Vel it 

[2] +01 IF L#0 

bea +0,2Z<R+1 

[L4J} £4:%0,2Z<(L-1)ACKER 1 

[5] £2:Z<(L-1)ACKER L ACKER R-1 
V 


Example : 


3 ACKER 1 
i 


10.6. Optical Illusions 


Some well known optical illusions are easy to effect on the screen 
of a small micro. Here are three of the best known simple illu- 
sions - realising them on a computer gives the opportunity to 
experiment with variations on a basic theme. The defined func- 
tions used are described in Appendix 1. 
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Which small square has the larger area? 


O° 0: 511 259 FLLLBOX 3: 5 pl00 100. 200.200 4 
17S A250: 7S: AIS: 
320 1259 375 L73* 1 


O 0 5141 255 BOXES 1 49 300 100 400 200 


Which horizontal line is the longer? 


A<2 REFLECT 1 REFLECT 3 3p0 0 
L 2 
di 
B<2 REFLECT 1 REFLECT 3 300 0 
1 2 
Ls 


FO OF OO © 


AL;2]<A[L;2]—-4 
BL;2]<+BL;2]+4 
“10 5 10 5 SKETCH A,[1]B 


Which circle is the largest? 


L3: 30> 0: 35-90 
t- A0: 50 
a 0 

Ie3: sip 30. <S0> 0 
50 50 10 
TO> SO 10 

W<0 0 120 120 


Use W SKETCH T & W CIRCLES Y then clear the screen and 
try W CIRCLES Y alone. 


Appendix 1 


Graphics 


There are many points in the preceding chapters, particularly the 
one in geometry. where the exposition of mathematics 1s greatly 
enhanced by a display of the results of running a program. 


APL itself contains no graphics primitives and so the production 
of good pictures depends on the availability of some additional 
resource — normally a separate graphics package — whereby 
arrays can be rendered into pictures. 


It is impossible to anticipate the range of graphics packages 
which might be used in conjunction with APL in the classroom, 
so it has been assumed that a basic package is available con- 
taining a number of primitive graphic functions which are 
described below. 


A common underlying concept is that of a “window.” Conceptu- 
ally there are two sets of co-ordinates which the user must keep 
in view, graphics co-ordinates and problem-space co-ordinates. 
The former are determined by the physical nature of the equip- 
ment and are intimated in the appropriate manual. To be spe- 
cific, suppose that the co-ordinates of a screen boundary are 
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0,199 319,199 


0,0 319,0 


The problem-space co-ordinates on the other hand are deter- 
mined entirely by the user, who perceives the screen boundary as 
determined by the extreme values of x and y which he/she 
wishes to be represented on the screen. 


X,Y 


X,Y 


If the user co-ordinates of the bottom left and top right corners 
of the screen are (x,,y,) and (x;,y>) then the vector 


X15V15Xo,Y2 
is called the “window.” 


A fundamental auxiliary function is XFM_ which transforms a 
matrix of values in problem-space co-ordinates into the corre- 
sponding values in graphics co-ordinates. The left argument of 
XFM is the window, the right argument consists of a single 
column of x-co-ordinates, followed by one or more columns 
which are interpreted as y-values. Given the screen dimensions 
above, XFM is 
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v Z<+L XFM R:I 
Ee T+«1 
[2] Z+((1+pR),1)0INT 319x(RE31)]-L201])+ 

—/L(3 1] 

L3] 21:70 IF ( 1+pR)<I<I+1 
C4] Z+Z,INT 199x(RCL3;TIJ-L0L2))+-/264 2] 
[5] -+L1 

V 


The function INT rounds a number to the nearest integer (see 
Appendix 2). The values 319 and 199 are screen boundary co- 
ordinates (see previous page) and should be replaced by appro- 
priate values for the user’s own system. The graphics functions 
themselves are 


IL DRAW fF 


L : window 
R : (NxP) data matrix, col | = x-values 
other cols. = sets of y-values 
Draws (P-1) sets of (N-1) connected line segments - the first 
given by ARL;1 2], the next by AL;1 3] and so on. 


i, DOT a 


As for DRAW except that the points are plotted but not con- 
nected. 


IL PIN B 
L : window 


R : 2 element vector (x,y) 
Draws the single point R as a pin-point or pixel. 


156 APL Programs for the Mathematics Classroom 


LZ JOIN R 


L : window 

R : 2 element vector (x,y) 
Plots a point and draws a line to it from the previously plotted 
position. 


IL SKETCH R 


L : window 
R : 3 column matrix. 
cols 2,3 are sets of (x,y) co-ordinates 
respectively; 
elements of col | are 0 = move, 1 = draw. 
Similar to DRAW for a single graph but adding a 
“penup/pendown” facility. 


L LINES FR 


L : window 
R : (Nx4) matrix. Each row 1s a set of 4 
co-ordinates in the form (X,,Y;,X>,Y>) 
representing the start and end points of 
a line. 
Draws as many lines as there are rows of R. 


LI BOXES & 


L : window 
R : (Nx4) matrix. Each row 1s a set of 
4 co-ordinates in the form (X,,Y},X>.Y>) 
representing the lower left and upper right 
corners of a box. 
Outlines as many rectangles as there are rows of R. 
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I FILLBOX Rk 


L : window 
R : (Nx5) matrix. First 4 columns as for 
BOXES. 5th. column is an integer 
representing the colour with which 
the box is to be filled. 0=no fill. 
Fills as many rectangles as there are rows of R. 


L ARCS fF 


L : window 
R : (Nx5) matrix. Each row is a set of 
4 co-ordinates in the form (X,,Y;.X2,Y>2) 
representing the end-points of an arc, 
followed by the angle in degrees 
of the arc connecting them. 
Draws as many arcs as are defined by the rows of R. 


Gh -CLAC LBS) At 


L : window 
R : (Nx3) matrix. Each row is (x,y,r) where 
(x.y) are the co-ordinates of the 
centre and r is the radius. 
Draws as many circles as are defined by the rows of R. 


HIST fh 
R : vector 
Draws a histogram the heights of whose boxes are the elements 


of R. The histogram is automatically scaled so that the picture 
fills the screen. 
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In addition it is useful to have a function REFLECT which for a 
matrix R, negates the values in the column(s) L, and joins the 
result to the original matrix 


REFLECT: Z<k & 
Zl ho eG ee ZR LZ 


Example : 


Ae2 ZO 46.56 
1 REFLECT fF 2 REFLECT F 


nfm £ 


3 
5 
“3 
es 


This set of functions covers all the graphics needs of the algo- 
rithms described in this book. The user must code them using the 
system functions available in the APL/graphics system to hand. 
perhaps adapting them to a set of functions which fit more 
naturally within that system. 


It should be added that if no graphics system is available, simple 
graphics may nevertheless be produced on a character-box basis 
as in the function LIFE (line 8) in Chapter 10. For example, the 
sequence 
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X16 
' *'£1+(0 24+110)°.=8' (X-2)xxX-4'] 


draws a graph of the parabola y = (x-2)(x-4) from x = | to 6. 


Appendix 2 


Idioms and Utilities 


This section contains a collection of general purpose programs 
relevant to mathematics. Some are useful in their “stand-alone” 
role - others occur as auxiliary functions in programs throughout 
this book. Many of the phrases given here have become part of 
the folk-lore of APL, and are recognised almost like primitives by 
habitual users. 


A2.1. Rounding, Averaging, and Removing 
Duplicates 
Round number(s) R to L decimal places 
ROUND : (10*-Z)xL.5+Rx10*Z 
L : integer (possibly 0 or negative - if negative 


R is rounded to nearest 10*—Z 
R : numeric scalar, vector or array 
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INT. eo 


Average of a vector R 


AVERAGE : (+/R)+pR 


Remove duplicate elements from a vector 
REMDUP : ((fif)=10R)/R 
alternatively 


REMDUP : (V/<\Ro.=R)/RB 


Remove duplicate rows from a matrix 
REMDUPM: (1 1&<\RA.=QR)4R 


A2.2. Sorting and Ranking 


Sort a vector R in upward/downward sequence 


SORTUP : RLA4R] SORTDWN : RLYER] 


Obtain upward/downward ranks for vector R 


RANKUP : AAR RANKDWN : AYR 


Appendix 2. Idioms and Utilities 163 


Obtain upward tied ranks for vector R 


TRANKU : .5x(AAR)+¥VAOR 


Obtain downward tied ranks for vector R 
TRANKD : .5x(AVR)+VVOR 


Obtain downward rank for vector R with equal members being 
given lowest possible rank (schoolmaster’s rank) : 


SRANK : (A4YR) (RIB) 
Change Y to 4 _ for upward schoolmaster’s rank. 


Examples : 


MARKS 
Of Gan 22 27: OF 22 
TRANKU MARKS 
he 3G. ain S. Ee. He ae 
TRANKD MARKS 
Sod! Dard. S495 se 
SRANK MARKS 
25d. 2 2. 


Merge vector R into vector L with items in correct order 
MERGE : (L,R)CAL,R] 
Example : 


0 10 20 30 MERGE 27 13 32 
0 10 13 20 27 30 32 
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Sort each row of a matrix in ascending order 
RSORTU : To(RLARI)CAl(AR<,R)+ 1+T«pR] 
Sort each row of a matrix in descending order 


RSORTD : To(RCYRI])CAT(YR+,R)+ 1+7+poR] 


Example : 

M 

ce oF ae 

2, & 

b. “8 “a 
RSORTU M 

7 10 

a « 2 


Sort (a) rows (b) columns of a matrix in ascending order 
UPMATR : RLAI(I/,2) 1883 J 
UPMATC : RL34l(1/,2) LF] 


n.b. leftmost column/topmost row is most significant within row/column 


Example : 
UPMATR M 
2. «Gi 35 
GS 3: 2 
oF £0 


This technique can be extended to sort the rows of a character 
matrix in alphabetic order : 
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UPSORT : UAVLUPMATR UAE] 


Example : 


MC 
SANTA 
CLAUS 
ESQ 
UPSORT MC 
CLAUS 
B'S@ 
SANTA 


A2.3. Statement Joining 


“Glue” 


In the introduction the use of the ‘&’ symbol was described as a 
means of keeping two related statements together on the same 
line. In some implementations this 1s acceptable as a correct APL 
construct; in others one must either write two separate statement 
lines, or else use a device such as the following to “glue” two 
statements together 


Pas 08 ae 


This is a technique which is sometimes frowned on by APL 
purists, because of a trap which can arise from e.g. 


A<2,00B<2 


Although apparently identical, B is a scalar whereas A is a 
l-element vector. Provided that this possibility is appreciated, 
the device can be used to compress two simple statements onto 
one line, and can be useful when the rightmost statement is an 
output message, e.g. 
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A<+2,0oLK'PROCESS TERMINATES ' 


Merging Branch Statements 


Another statement joining technique uses the fact that a branch 
to a non-null vector means branch to its first item, ignoring all 
other items. Thus 


>hi.Z< ... 


joins a branching statement to another statement on its right. 
The result of the other statement must be a numeric scalar or 
vector - sometimes it may require an extra ravel to achieve this, 
1.€. 


SO le ie Gx 


Read this as “branch to L1, having assigned Z to ....” See SER] 
in Chapter 4 for an example of this technique. 


Mid-line Assignment 


Yet another device which is frequently used to merge two or 
more statements on one line is mid-line assignment, of which the 
technique described above is a special case. It is frequently used 
to increment loop-counting variables (see line 2 of AMICNOS), 
and another general example from Chapter 2 is 


PRIMES : (2=+/0=Te°.|7)/T«.ik 
It is a good idea to read this as 
“C ....) reduce 1 which is incidentally assigned to T.” 
Some APL users discourage mid-line assignment on the grounds 
that in reading the line above from left to rnght, you cannot 


“understand” T in the parentheses until you get to a later stage 
in the APL sentence. While it is true that overuse can lead to 
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incomprehensible program lines, its use in moderation can often 
be justified. 


A2.4. Branching and Prompting 


Two other useful functions are IF and ASK. The function IF 1s 
used repeatedly in this text. ASK is used occasionally to show its 
possibilities; in classroom practice it can provide user-friendly 
interaction with some of the mathematical functions. In order to 
keep the mathematical aspects as clear as possible, such “cos- 
metics” have been trimmed in the programs given. This 1s not 
however to underestimate their importance in the overall presen- 
tation of mathematical software. 


IF : R/L 


L : label 
R : condition 


This construction is widely used in this text to clarify conditional 
branches by writing e.g. 


+01 IF R<I<I+1 rather than >(f<7<7+1)/51. 


Prompting 
V Z<ASK Ff 
[1] R 


FD 5 >D1 IF 0#pZ+{J 
[3] +0,2<10 
[4] £1:2<#Z 

Vv 


Example : 


A<,ASK 'GIVE 2 NUMBERS: ! 
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following which A will have the value of whatever APL 
expression (which includes of course a constant scalar or vector) 


is given by the user in response to the prompt. See CAST9 in 
Chapter 2 for a further example. 


A2.5. Matrix Manipulation 


Column and row addition of vector L to matrix R 
COLADD : R+(pR)oL 
ROWADD : &((%pR) pL)+R 
L : vector with same number of items as R has 
cols/rows 
R : numeric matrix 


Example : 


O 1 2 ROWADD 3 3019 


1, 2 2 
2 o 
o. AO 1S: 


Upper and lower triangular RxR matrices 


These have 1’s in upper/lower triangle (excluding the leading 
diagonal), 0’s elsewhere. 


UPT : (i1R)e°.<iR and LOT : (1f)°.>1iRf 
Ditto including the leading diagonal 


UPTD : (1f)°.<if and LOTD : (1Rf)°.2if 
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OOR 
Orb 
bp BE 


Unit matrix of order R 


UNIT :(1f)°.=18 or UNIT : (RVR) oCR+1)t1 


Border a numeric matrix with zeros 
BORDER : (-—2+pRf)+(1+oR) +R 


A2.6. Replication 


The following function implements the function REPLICATE 
which is to be found in some APL’s as an extension of COM- 
PRESS. R is a vector (character or numeric), and L is a numeric 
vector of the same length indicating for each item of R how 
many times (possibly zero) it is to be replicated. COMPRESS is 
then a special case in which only 0’s and 1’s are allowed in L. 


v Z<«L REPLIC R 
[1] Z<+10 
[2] 21:70 IF 0=pL 
[3] Z<Z,(1tL)pitR 
C4] I<ivL & R<1ivR 


So sed 
Vv 
Example : 


23 0 1 REPLIC ‘ABCD ' 
AABBBD 
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A2.7. Without 


Another function which is to be found in some APLs as dyadic 
“not” is WITHOUT, which is approximately modelled by 


WITHOUT : (~LeR)/L 
Example : 


(1140)WITHOUT 4 5 6 
123 7:8 9 10 


A2.8. Bit Manipulation 


In the following set of idioms involving the scan operator, a bit is 
considered to be switched “on” if it has the value 1, “off” if it 
has the value 0. 


Switch on all bits from and after the first | : V\R 
Switch on the first 1, off all bits thereafter a» SNA 
Switch on all bits after the first 0 . S\f 
Switch off all bits after the first 0 > A\R 
Count the number of 1’s before the first 0 . LANE 
Examples : 

M1 


PER OOO 
LF OORKO 


FO FPF OO F © 
Fe OO FP FP © 
oO OF CO FF 
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V\M1 A\M1 
0 0 0 al 000 0 0 
O1t114%1t1 0000 0 
Ord, ds eb oh 000 0 0 
111411 100 0 0 
i? ae, el A Ss 1000 0 
A, A he he 113110 

<\M1 <\M1 
OOF 20 10 3 O11%11 
O10 0 0 O1i14%11 
O10 0 0 Oe Ae od, Ay 
A. Q.0:. Or <0 TO 2. aL 
A. 20) Or iO OD dO A A 
1000 0 ais Sate fale “aes 


A2.9. Some String Handling Functions 


Although mathematics makes less demands than most other 
application areas on the string handling capabilities of APL, 
these are occasionally useful, and so three general utilities are 
given here. 


Blank out all occurrences of a specified set of characters 


BLANK : T\(T<~LeR)/L 


L : string in which replacement is to take place 
R : string of characters to be replaced by blanks 


Example : 


'(2X-3Y).(4X+5Y)' BLANK 'XY().' 
2 3 ae 
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Find all occurrences of one string in another 
FINDALL: RA.=( 1+1pR)$((pR), pL) oL 


L : string to search in 
R : string to search for 


The result of FINDALL is a bit vector of the same length as L, 
with 1’s marking the start position of matches of R within L. 


Example : 


'(2X-3Y).(4X+5Y).(X-7Y)' FINDALL ').(' 
OP Oe Ov 0). “Qe 2.2 0e-0. 20-0. 0 0 <0" 0 0) O-0: 20: 0 


A more general utility which covers many of the circumstances 
where string handling is useful in a mathematical context 1s 
REPLACE which replaces all occurrences of one string.in L with 
another. 


V Z2<L REPLACE A37T371;723U 
[1] T1<(7<Ri's')VR & T2<(T-1)+h 
[2] Z<e''  & U<L FINDALL T2 
[3] 61:22 IF 0=v/U 
[4] Z<Z,( (T+ 14Ui1)4L),71 
Roa L<(T<T+e7T2)+L & U<TYU 


[6] >D1 
[7] £2:2<Z,L 
V 


L : string in which replacements are to be made 
R : string to be replaced, 

followed by semi-colon, 

followed by replacing string 
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Line 1 separates out the replaced and replacing strings; line 2 
establishes the points where replacement is to take place, and the 
rest of the function is a loop to do the replacements one by one. 
If ’;’ is one of the characters which has to be replaced, a different 
separator character must be used, and an appropriate adjustment 
made to line 1. 


Examples : 


'2 -3 4 +5!' REPLACE '-; ' 


2 3 4 +5 

'2,399,999' REPLACE '99;00' 
2,300,009 

'2,399,999' REPLACE ',;3' 
2399999 


A2.10. Testing for Numeric/Character 


The function TESTNUM returns | if R is numeric and 0 if R 1s 
character 


TESTNUM : O=1t0pR 


As noted earlier, the programs in this collection have been 
stripped of data-checks in order to help clarify essential program 
structure. For versions intended for routine use, and which 
might be passed on to other users, it is wise to commence each 
program with a prologue which does routine checking. For 
example, if both arguments of a function must be numeric, L 
must be a scalar and R a matrix, begin the function 


[1] +L0 IF (TESTNUM L)A(TESTNUM R)A 
(O=ppL)A2=ppF 

[2]. 405, 0ple'2RROR' 

[3] ZO: ... function proper 
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A2.11. Timing Function Execution 


One utility function which most users eventually want is a 
TIMER function. The exact contents depend on how the user’s 
own system reports CPU time used. Here is a TIMER which 
assumes that current CPU time is recorded in DATL2], and that 
R is the character string form of a phrase which it is required to 
time for L executions 


V Z2<L TIMER B31 
[1] T<+0 & ZHIATL2] 
[2] £1:302 IF L<I<I+1 


[3] oR 

[6] >D1 

[7] £2:2Z<QATL2)]-Z 
V 


Example (on one specific computer) : 


100 TIMER 'X<?1000100' 
86 


Appendix 3 


Graphics Functions in |-APL 


Using I-APL and the supplied workspace PGRAPH, the fol- 
lowing functions implement the graphics functions of Appendix | 
on IBM PC’s and compatible machines. Direct definitions should 


be entered as given, 1.e. a and w correspond to L and R else- 
where in the book. 


DRAW: (aDRAW O 14w), 
(32767 32767 0 1-0-0 0-0, 1+41%0@) 
PGPLOT LaXFMw[3:1,7] : 1=T+ 1tow :10 


PLOT: (aPLOT 0 14¥w), 
PGPLOT LaXFMw[l31.7] : 1=T+ 1tow :10 


PIN: PGPLOTaXFM1 2peuw 


V L SKETCH R;T 
[1] 01:70 IF O021+0R 
Ea Al13;1J<1 
[3] T+ 1+RO31)10 
C4] ZL DRAW RliT;2 3] 
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[5] R+(T,0)VR 
[6] D1 
V 


The PGRAPH functions have no concept of “current graphics 
pointer,” and so the function JOIN must be programmed sepa- 
rately for each application, with a variable recording the co- 
ordinates of the last plotted point. The functions for LINES, 
BOXES, and FILLBOX are 


EINES :32767 32767 0 1-00. “4. “4: toa 
PGPLOT LaXFM(2 .5xow) pw 


v L BOXES R 
[1] £1:70 IF O021+pR 
[2] E DRAW AIS se 3-4 3 aS 
eee ye Ue Oe 
[3] R<1 OVE 
C4] +D1 
V 


v GL FILLBOX &;T;U 
[1] £1:70 IF 021+ oR 
En T+LL XFM 2 2oR[13] 
[3] U+|-/TL32] 
C4] TT Oe0) eT {129 ] 
Basi (32767 32767 0 1 0,U,0 0 1)PGPLOT T 
[6] >D1 
V 


The function to implement circles assumes that each circle is 
drawn as an 18-sided polygon. If higher resolution is required, 
change the ‘19’s’ in line 2 to one more than the number of sides 
in the approximating polygon. 
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VY i; CIRCLES fF 
[1] 21:70 IF 021+poR 
[2] L DRAW (19 29R[L1;12]+ 
Rl13;3] CIRC 20x19 
[3] R+1 OVR 
ES +D1 
V 


For CIRC see Section 6.1. 


ARCS defines arcs as a series of connected line segments. The 
function ARC embodies the necessary co-ordinate geometry to 
give the anti-clockwise arc of angle R[5] degrees connecting the 
points RL1 2] and RL3 4]. To obtain the clockwise arc reverse 
the pair of co-ordinates. 


The function ARC assumes that points will be connected at 5 
degree intervals. If greater or lesser resolution is required, change 
the value 5 in line 2 of ARCS to the desired value. 


The auxiliary function PTS requires for R an anti-clockwise 
range of degrees as a 2-item vector, and for L the resolution in 
degrees as described above. 


V L ARCS R 
[1] 21:70 IF O21+pR 
[2] ZL DRAW 5 ARC RL[1:3] 
[3] +1 OVP 
C4] +D1 
V 
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ee 
[2] 
[3] 
[4] 
ES 
[6] 
[7] 
[38] 
[9] 
Laos 


a 
[2] 
[3] 
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VY 2<L ARC R3T3U3V3W3X37Y 
U+.5x(O 1+R)+180 
T<—f#Y<2 2044+R 
W+((+/7*2)*.5)#2x10U 
>f1 IF 04147 

V+o.5 

+L2,00X<—-x1+T 


[1:V+ 30:+/oT7 


X*+—x14*7 


L2:Z«(AV Y)+(XxWx20U)x1 “1x1 20V 


Z+((o7)0Z)+7«W CIRC L PTS(U-2x 14*T), 
U+(180+Xx90)++/7<(180+01)xV,U 


V Z2<PYS PR 

R<«R+0,360x>/R 
Z+Lx(LAC1)]+L2)+1L(-/OR) +2 
Z+360|R(1)],Z,R02)]xiRL2Je% 142 
Vv 


AV: (+4w) #1t ow 


The following function implements HIST in Appendix 1: 


HIST:(0 0 0 1,7,32767, (T«L600+1+pw), 
0 0,0w)PGPLOT( (pw) ,1) 9 Lwx180+1/w 
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Ackerman’s function, 150 

Adaptive Simpson’s rule, 90, 
129-130 

Addition of polynomials, 32 

Amicable numbers, 14 

Angle between lines, 82 

Area of a triangle, 59 

Arabic numbers, 18 

Area of circle, 55 

Argument of vector, 81 

Arithmetic progressions, 33 

Astroid, 71 

Average, 103,162 

Axis of graph, 8 

Backward counting, 23 

Best fitting poly, 113 

Beta function, 126 

Binary fractions, 17 

Binomial coefficients, 43 

Binomial distribution, 99 

Birthday problem, 102 

Bisection method, 91 

Bit manipulation, 170 

Bordering matrix, 169 


Box-Muller formula, 106 
Branching, 167 

Buffon’s needle, 108 
Cardioid, 75,76 

Cartesian co-ordinates, 82ff. 
Casting out nines, 21 
Centigrade degrees, 55 
Chessboard, 146 
Chi-squared distribution, 129 
Chi-squared test, 125 
Circle, 67 

Coins, tossing, 108 
Combinations, 140 
Complex numbers, 24 
Complex roots of unity, 25 
Compositions, 142 
Compound interest, 56 
Computer arithmetic, 22 
Conditional branching, 167 
Conic sections, 67,72 
Continued fractions, 51 
Co-ordinate geometry, 81 
Convergence, 42 
Correlation, 11S5ff. 
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Cos formula, 60 
Covariance, 115ff 

Cube, 79 

De Moivre’s theorem, 35 
Decoding, 19 

Deltoid, 76 

Derangements, 138 
Descriptive statistics, 103ff. 
Determinant, 31 

Dice, throwing, 108 
Difference of squares, 27 
Differences of series, 45 
Differential equations, 96 
Differentiation, 32 

Digit sum, 21 

Direct definition, 4 
Duplicate removal, 162 

E, series for, 50 

Ellipse, 67,73 

Encoding, 19 

Envelopes, 71 

Epicycloid, 68,75 

Euclid’s algorithm, 15 
Euclidean norm, 62 
Euler’s series, 47 

Euler’s method, 96 

Even integers, 7 
Expansion of brackets, 28 
Exponential distribution, 106 
F distribution, 130 
Factors, 13 

Fahrenheit degrees, 55 
Fibonacci series, 40, 46 
Finding substrings, 172 
Formulae, 55ff. 

Forward counting, 23 
Fractional parts, 19 
Frequency distribution, 109 
Gamma distribution, 129 
Geometric progressions, 33 
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Glue, 165 

Goodness of Fit, 125 

Graphics co-ordinates, 153 

Gregory’s series, 47 

Hastings formula, 127 

Highest common factor, [5 

Histogram, 157 

Hyperbola, 67,73 

Hypergeometric distribution, 
101 

Hypocycloid, 68,75 

I-APL Graphics Fns, 175 

Idioms, 161 ff. 

Images, 81 

Integration, 32,87 

Interpolation, 52 

Intersection, 34 

Inverse of matrix, 30 

Isomorphisms, 11 

Iteration method, 92 

Joining statements, 165 

Kendall coefficient, 119 

Koehler’s formula, 128 

L-norm, 63 

Lagrange’s formula, 52 

Latin squares, 143 

Least squares fitting, 113 

Lexical ordering, 135 

Life, game of, 148 

Limacon, 76 

Logarithm tables, 11 

Longest journey, 61 

Lower triangular matrix, 168 

Lowest common multiple, 16 

Magic squares, 144 

Mann-Whitney test, 123 

Mastermind, 147 

Matrix inverse, 30 

Matrix manipulation, 168 
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Matrix operations, 30 
Matrix product, 30 
Matrix quotient, 30 
Median, 104 
Mensuration, 20 
Merging vectors, 163 
Mid-line assignment, 166 
Mid-point method, 97 
Mode, 105 

Money conversion, 20 
Monotonicity, 42 
Moran’s formula, 127 
Mortgage repayment, 58 


MultipIn of polynomials, 32 


Multiplication tables, 10 


Negative exponential dist., 106 


Nephoid, 75 
Neville’s algorithm, 53 


Newton-Raphson method, 


94ff. 


Non-parametric tests, 117ff. 


Norms, 62ff. 


Normal distribution, 126ff. 


Number bases, 16 
Numerical integration, 87 
Odd integers, 7 

Optical illusions, 150 
Parabola, 68,74 
Parametric equations, 67 
Partition values, 104ff. 
Partitions, 143 

Pascal’s triangle, 44 
Patterns of points, 84ff. 
Pen-up/pen-down, 66 
Percentiles, 104 

Perfect numbers, 14 
Permutations, 135ff. 
Perpendicular distance, 81 
Perspective drawing, 80 


Pi, 47-49 

Poisson distribution, 100 

Polar co-ordinates, 82ff. 

Polynomials, 31,36 

Polynomial interpolation, 
53 

Positive integers, 7 

Precision, 23 

Present value, 58 

Primes, 12 

Probability density, 126ff. 


Problem-space co-ordinates, 
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Product moment correlation, 


115 
Prompting, 167 
Pythagoras theorem, 62 
Pythagorean triples, 63 
Quadratic equations, 29 
Random numbers, 106ff. 
Range, 105 
Rank correlation, 118 
Ranking, 162ff. 


Recurrence relations, 39ff. 


Recurring decimals, 16 
Recursion, 5,46,149 
Regression, 113 
Removing duplicates, 162 
Replacing strings, 172 
Replication, 169 
Residuals, 114 

Roman numbers, 18 
Root finding, 91 ff. 

Root 2, series for, 50 
Roots of polynomials, 36 
Roots of unity, 25 
Rotation, 77ff. 
Rounding, 161 

Runs test, 117 
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Sample sizes, 131 

Scaling axes, 8 

Scatterplot, 112 
Schoolmaster’s rank, 163 
Selections, 141 

Series, 39ff. 

Sets, 34 

Shortest journey, 61 

Sign test, 120 

Simpson’s rule, 89 
Simulation, 1O8ff. 
Simultaneous equations, 30 
Sin formula, 60 

Sliding ladder curve, 71 
Sorting, 162ff. 

Spearman’s coefficient, 118 
Square numbers, 8 
Standard deviation, 104 
Statistical tables, 126ff. 
Stem and Leaf Plot, 110 
String handling, 171ff. 
Student’s t distribution, 128 
Subsets, 35 

Sum of squares, 28 
Supersets, 35 

Surface area of cone, 56 
Switching bits, 170 

T distribution, 128 

Table lookup, 10 

Taxi-cab metric, 63 


Test for numeric, 173 

Tied rank, 163 

Timing function, 174 

Tower of Hanoi, 149 

Transformations, 77 

Trapezium method, 97 

Trapezium rule, 89 

Triangle formulae, 59 

Triangular matrices, 168 

Triangular numbers, 8 

Trigonometric functions, 51 

Trigonometric tables, 11 

Two-way frequency, I11 

Uniform distribution, 106 

Union, 34 

Unit matrix, 169 

Upper triangular matrix, 168 

U-test, 123 

Utilities, 161 ff. 

Variance, 103 

Velocity, 56 

Volume of cylinder, 56 

Volume of sphere, 55 

Whist, 145 

Wilcoxon Test, 122 

Wilson-Hilferty approximation, 
129 

Window, 153 

W-test, 122 

Yates’ correction, 125 
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Index of Programs 
and Variables 


ACKER 150 BINPERM 138 COMPOSN 143 
ADSIM 91 BINPROB 99 COMPS 143 
ALLS 35,36 BISECT 92 CONFRAC 51 
AMICNOS 15 BLANK. 172 CONIC 72 
ANG 82 BORDER 169 CORM 115 
ANGLEA 59 BOXES 150,156,176 COS 51 
ANTILOG 11 BTONUM 140 COVM 115 
AP 33 BTONUMI 142 CPCT 129 
ARABIC 18 BUFFON 108 CPDF 129 
ARC 15], 178 BWD 24 CRU 25 
ARCS 157,177 CAST 22 CTOP 82 
AREAS 88 CAST9 22 CTOP!1 83 
ARG 81] CHBOARD 146 CUBE 79 
ASK 167 CHVECT 146 DECFRAC 18 
ASTROID 71 CCOF 115 DERANGE 138 
AV 103,178 CHISQ 125 DET 31 
AVERAGE 162 CI 56 DICE 108 
AXISLAB 8,87 CIRC 67 DIF 45 

BASE 17 CIRCLES 66,157,177 DIFFSQ 27 
BETA 126 COINS 20 DIGSUM 21 
BINCOEF 43 COLADD 168 DIST 81 


BINFRAC 17 COMBS 140 DIV 24 
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DRAW 65,155,175 


El 50 

E2 50 
EHPARM 69 
EHCURVE 70 
ELLIPSE 67 
EPIHYP 75 
EST 113 
EUC 15 
EUL 96 
EULER 47 
EVAL 17 
EVENFAC 14 
EXPAND 28 
EXPEC 125 
EXPRESS 17 
EXTEND 76 
FACS 13 
FACTINT 135 
FACTORS 13 
FIB 46 

FIB1 46 
FIB2 47 


IF 167 
IMAGE 81] 
INRANGE 112 
INT 162 
ITEROOT 93 
IX 34 

JOIN 66,156 
KENDALL 119 
LAGR 53 
LATIN 144 
LCM 16 
LEXPERM 137 
LIFE 148 
LIMACON 77 
LINES 66,156,176 
LNORM 62 
LOG 11 
LONGEST 61 
LOOKUP 10 
LOT 168 
LOTD 169 
MADJ 103 
MAGIC 144 


FILLBOX 150,157,176 MATCHES 102 


FINDALL 172 
FPDF 130 
FREQ 109 
FWD 23 
GOFIT 125 
GP 33 

GPDF 129 
GREGORY 47 
HAND 145 
HANOI 149 
HATS 139 
HCF 15 

HIST 157,178 
HYPERB 67 
HYPGEOM 101 


MATCHPR 102 
MAXI 109 
MDICE 108 
MEDIAN 104 
MERGE 163 
METRIC 20 
MIDPT 97 
MM 147 
MODE 105 
MONOD 42 
MONOI 42 
MORT 58 
MUL 24 
MULTAB 9 
NEV 54 
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NEWRAP 94 
NEWRAPI1 95 
NINT 127 
NORM 62 
NPCT 127 
NPDF 126 
NUMDER 139 
ODDS 7 
OFF 54 
PADD 32 
PARAB 68 
PARENV 74 
PAS 44 
PASN 44 
PATTERN 84 
PCOFS 36 
PCTILE 105 
PDIFF 32 
PDRAW 114 
PERFNOS 14 
PERM 137 
PERMS 137 
PERSPEC 80 
PIN 155,175 
PINTEG 32 
PIl 48 

PI2 49 

PI3 49 

PI4 49 

PIS 49 

PLOT 65,155,175 
PMULT 32 
POISSON 100 
POLI 54 
POLY 32 
POLYFIT 113 
POWERS 9 
PRECTST 23 
PRIMES 12 
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